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Many biological, geophysical and technological systems involve the transport of 
resource over a network. In this paper we present an algorithm for calculating the 
exact concentration of resource at any point in space or time, given that the re- 
source in the network is lost or delivered out of the network at a given rate, while 
being subject to advection and diffusion. We consider the implications of advec- 
tion, diffusion and delivery for simple models of glucose delivery through a vascu- 
lar network, and conclude that in certain circumstances, increasing the volume of 
blood and the number of glucose transporters can actually decrease the total rate 
of glucose delivery. We also consider the case of empirically determined fungal 
networks, and analyze the distribution of resource that emerges as such networks 
grow over time. Fungal growth involves the expansion of fluid filled vessels, which 
necessarily involves the movement of fluid. In three empirically determined fungal 
networks we found that the minimum currents consistent with the observed growth 
would effectively transport resource throughout the network over the time-scale 
of growth. This suggests that in foraging fungi, the active transport mechanisms 
observed in the growing tips may not be required for long range transport. 

Keywords: Transport networks; fungal networks; vascular networks; advection- 
diffusion. 
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1 Introduction 



Many biological, geophysical and technological systems involve the transport of 
material over a network by advection and diffusion ifTTl fT4l l38l l42l l52l . and it 
is common that this material can leave, decay, be lost, consumed or delivered 
as it propagates. Indeed, fluid transport systems are found in the vast majority 
of multicellular organisms, as the component cells of such organisms require re- 
sources for metabolism and growth, and diffusion alone is only an effective means 
of exchange at microscopic length scales BUI . Molecules of interest are car- 
ried by advection and diffusion through the cardio-vascular networks of animals 
E@l22l3ai2ai23mi2ai§a,the mycelial networks of fungi flSHHl, the xylem 
and phloem elements of tracheophytes (vascular plants) Il44l I5T1 |59l , and vari- 
ous body cavities of many different animals. For example, oxygen is transported 
through the lungs of mammals and the trachea of insects, while suspension feeding 
animals (including sponges, clams, brachiopods, many arthropods, fish, ascidians 
and baleen whales) pass water through various chambers of their bodies, captur- 
ing the organic particles that are needed for survival [40]. Similar mechanisms 
of transport are also found in geological and technological systems, such as rivers 
and drainage networks [6], gas pipelines, sewer systems and ventilation systems 

nana. 

In all of these cases the particles of interest diffuse within a moving fluid, which is 
constrained to flow within a given network. The bulk movement of fluid is referred 
to as advection, convection or mass flow, and in general the fluid in question travels 
with a mean velocity that varies over the network. The mean velocity of fluid flow 
may vary by several orders of magnitude, as, for example, the velocity of human 
blood drops from lm s _1 in the aorta to around 1mm s _1 in the capillaries lfT5ll2TTl . 
Given a network and a distribution of velocities, we may wish to calculate how an 
initial distribution of resource changes over time. For example, we might want to 
know how a patch of pollutant will spread within a drainage network IfTTl 14211521 . 
how a drug will spread within the cardio-vascular system Il8ll^ l2l^l30ll53ll54ll60ll . 
or how nutrients will be translocated within a fungal network |[T6l [36l . In this 
paper, we consider the particular cases of modelling the delivery of glucose via a 
vascular network, and modelling the translocation of nutrients in a fungal network. 

Koplik et. al. [39] describe an effective method for calculating the exact moments 
of the transit times for a neutral tracer across an arbitrary network that contains 
a flowing medium, but which initially contains no tracer. We have advanced their 
methods to handle resources that may be consumed or delivered out of the network, 
while the resource that remains in the network moves by advection and diffusion. 
More specifically, we suppose that each edge in the network has a local deliv- 
ery rate Rij , which represents the probability per unit time that any given unit of 
resource will be consumed, lost or delivered out of the network. The effect of 
including a delivery term can be significant and somewhat counter-intuitive: we 
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will see, for example, that there are circumstances in which increasing the num- 
ber of blood vessels in a region can actually decrease the amount of glucose that 
is delivered to that region (Section [3]>. This problem is of particular bio-medical 
interest, as glucose delivery is essential to the survival of tumours and healthy tis- 
sue |[T7l[T9l[38ll54ll60l . As we shall see, to appreciate how the number of blood 
vessels in a region effects the total rate of glucose delivery, it is essential that we 
consider both the rate of delivery of resource out of the network and the topology 
of the transport network itself. 

To enable the assessment of the transport characteristics of arbitrary networks, with 
velocities that may vary over several orders of magnitude, we have developed a 
mathematical methodology that operates in Laplace space. We were initially mo- 
tivated to develop this algorithm by our interest in fungal networks. Peculiarly, 
the translocation of resource within fungal networks is much less well studied than 
transport in the other major multi-cellular kingdoms of life, but the ability of fungal 
colonies to translocate resources is ecologically critical ||6T1 . The relative roles of 
mass-flows (advection), diffusion and active transport are very poorly understood. 
Independent of exclusively fungal questions, fungal systems have the benefit that 
the network is accessible, and development can be readily followed through a se- 
quence of images. 

We have structured this paper to bring out the applications of our approach. As a 
consequence, a good part of the mathematical detail is in the Appendix. Although 
an important part of the paper's new results and machinery is in the Appendix, fa- 
miliarity with that material is not needed to understand the results that we discuss 
in the main text. The outline of this paper is as follows: Preliminary assump- 
tions and the fundamental equations governing advection, diffusion and delivery 



are discussed (Sections 2.1 and 2.2 1, and we stress the importance of the relevant 



time-scales for advection, diffusion and delivery (Section 2.3 1. We have developed 
a mixed method that enables us to calculate the exact concentration of resource at 
any point in space or time in an arbitrary network, and in the Appendix we de- 
scribe two efficient algorithms for updating the concentrations in a network over 
time, given an arbitrary, stepwise constant initial condition, and any number of 
point sources. In Section [5] we give a brief account of the convenience of solving 
the fundamental equations in Laplace space, and outline the key ideas and equa- 
tions of our approach. Alternative methods for solving the fundamental equations 
are outlined in Sectionf 



Finally, we apply our algorithms to a number of test cases, including a model of 
glucose transport in an idealized vascular network (Section [3]). We are motivated 
to understand how the geometry of a vascular network impacts upon the total rate 
of glucose delivery, as the effective use of anti-angiogenic drugs depends upon 
understanding the relationship between vascular pruning and nutrient delivery. We 
also apply our algorithm to a model of resource translocation across empirically 
determined, growing fungal networks (Section [4]). We note that changes in fungal 
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volume requires the movement of fluid: for example, the cytoplasm in a growing 
hyphal tube moves forward with the growing tip fill . In Section |4] we use our 
algorithm to investigate whether these growth induced currents are sufficient to 
supply the tips with the resources they require. In three empirically determined 
fungal networks we found that the minimum currents consistent with the observed 
growth would effectively transport resource from the inoculum to the growing tips 
over the time-scale of growth. This suggests that the active transport mechanisms 
observed in the growing tips of fungal networks may not be required for long range 
transport. 

2 Further details 

2.1 Preliminary assumptions 

We are interested in calculating the distribution of resource across a network of 
tubes, where the resource in question has a molecular diffusion coefficient D m , 
and where we are given four essential properties for each edge in the network (see 
Fig. 1). The edge connecting nodes i and j has: 

1. A cross-sectional area, denoted Sij(t). We assume that £y(i) is piece- wise 
constant in time, though in the Appendix we also consider the more complex 
case where Sij(t) varies continuously. 

2. A length, denoted 1^. As the location of the nodes does not vary over time, 
lij is constant. 

3. A mean velocity, denoted ity(i). This represents the mean velocity of the 
fluid in the edge, and we say that mj (t) is positive if and only if the current 
flows from node i to node j (so Uij(t) = —Uji(t)). By assumption, for each 
edge ij, Uij(t) is piece-wise constant in time. 

4. Finally, we suppose that resource in edge ij is delivered out of the network at 
a rate Rij, so if a particle is in ij for a short period of time At, the probability 
that it is delivered out of the network in that time is Rij At. 

While there is a single value for the molecular diffusion coefficient D rn , the disper- 
sion coefficient Dij(t) may be different for each edge. The value of Dij(t) reflects 
the tendency of adjacent particles to spread out within ij: they not only diffuse 
along the length of the transport vessels that comprise the edge ij, but also diffuse 
between the slow moving fluid by the edge of the vessels, and the relatively fast 
moving fluid in the centre of each vessel. 

If we consider the case where each edge ij is composed of some number of cylin- 
drical tubes of radius m (see Fig. [2]), and if the Reynold's number is small, we 
can calculate Dij(t) by using Taylor's dispersion coefficient for laminar flow in a 
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Figure 1 : Properties of a single edge in a resource distribution network. Sij (t) 

denotes the cross-sectional area of edge ij at time t, hj denotes the length of the 
edge, resource and medium flows along the edge with a mean velocity Uij{t), and 
resource is delivered out of the network at a rate Rij. Note that resource travels 
along each edge (and into other edges) by advection and diffusion, but the total 
rate at which resource in the edge is delivered out of the network is simply 
times the quantity of resource present in the edge. Also note that we do not need to 
assume that the edges in our network are straight, but we do assume that a single 
length scale lij captures the distance that particles must travel to move from i to j. 

cylindrical tube 11581 . This formula tells us that 

2 

Dij (t) = D m + Uij (t) 2 ^-. (1) 

In the case of a vascular network is simply the lumen radius of the edge ij, so 
we have Sij = 7rrf -. In plants, fungi or neural tissue each edge in the transport 
network can be modelled as a bundle of cylindrical tubes; in which case is the 
characteristic radius of the component transport vessels, and Sij is the total cross- 
sectional area of the transport vessels. 




Figure 2: Properties of an arbitrary resource distribution network. Each edge 
in the network is comprised of a single vessel or a bundle of transport vessels, and 
each edge has a length l^, a total cross-sectional area Sij(t), a mean velocity of 
flow Uij(t) and a local delivery rate Rij. Each edge also has a dispersion coeffi- 
cient Dij(t), as described by Equation (jl]). Note that the values of Dij(t) depend 
on the molecular diffusion coefficient D m , the velocities Uij(t) and the radius of 
the transport vessels within the edge ij. The nodes represent the point of contact 
between the edges: we assume that there is perfect mixing at each node, and we 
require a consistent concentration at node i whether we consider it to be one end 
of edge ij, or one end of any other edge connected to node i. 
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2.2 Fundamental equations 



We suppose that resource is lost or delivered out of the network at a given local rate, 
while the resource that remains within the network moves by advection and diffu- 
sion. Such a process will result in a spatial distribution of resource that changes 
over time. We only consider longitudinal coordinates along the edge ij, using 
real numbers x to denote distances from node i, where < x < kj. Each edge 
contains a quantity of resource, which must satisfy the one-dimensional advection- 
diffusion-delivery equation 

^ + ^ + ^^-A^ = 0, (2) 

where qij is the quantity of resource per unit length, Uij is the mean velocity, Dij 
is the dispersion coefficient and Rij is the rate at which a unit of resource is lost, or 
delivered out of the network. In other words, at time t and location x, the amount of 
resource in a Ax long slice of the edge is qij(x, t)Ax. The distribution of resource 
within each edge will vary over space and time, but if there is no direct link between 
the nodes i and j, we let Sij(t) = and q%j(x, t) = 0. This ensures that the sums 
in the following equations are properly defined for all pairs of nodes i and j. 

We wish to find the quantity of resource per unit length qij{x, t) at a given time t. 
A fundamental assumption underpinning the algorithms described in the Appendix 
is that there is perfect mixing at the nodes. In other words, the edge ij is only 
affected by the rest of the network via the concentrations at nodes i and j. This is 
only a reasonable assumption if the volume of the intersections between edges is 
negligible in comparison to the volume of the edges themselves. In effect, we as- 
sume that the nodes have an infinitesimal volume, so there can be no concentration 
gradients or boundary layer effects at the junctions between the edges. 

Crucially, the concentration at node i must be consistent across the edges ij, ik, 
etc, and we let Ci(t) denote its concentration at time t (amount per unit volume). In 
other words, for each edge ij we have 

where Sij(t) denotes the cross-sectional area at time t. 

It follows from our assumptions that the concentration profile in edge ij is com- 



pletely determined by Equation (23 1 together with the initial condition qij(x, 0) 
and the boundary conditions Sij(t)ci(t) and Sij(t)cj(t). By Fick's Law the rate at 
which resource leaves node i along edge ij is given by 



Jij(t) 



u ij (t)q ij (x,t) - Dij ^Q^' t 



(4) 

x=0 
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We assume that resource cannot accumulate at the nodes (as they have zero vol- 
ume), so any resource that enters node i along edge ij must leave node % along 
some other edge ik. Our framework can accommodate the case where resource is 
introduced at node i at some given rate Ii(t) > 0. However, if node i is not an inlet 
node (that is, a point where resource enters the network), we have Ii(t) = 0. In 
either case, Equation Q implies that the net rate at which resource leaves node i is 

hit) = 



Uij(t)qij(x,t) - D, 



dqij(x,t) 
' dx 



(5) 



Equations Q and ( 24 ) describe the current of resource, but we can also consider 



the current of fluid passing through a given point. Henceforth the term current is 
reserved for the quantity of resource that passes a given point per unit time, while 
medium-current refers to the volume of the advecting medium that passes a given 
point per unit time. The medium-current in edge ij is simply Uij(t)Sij(t), so the 
net medium-current leaving node i is 

Fi(t) =^2iHj(t)Sij(t). (6) 



2.3 Critical time-scales for advection, diffusion and delivery 

For an edge of length I, mean velocity u > 0, dispersion coefficient D and local 
delivery rate R, there are three critical time-scales: 

I 

tA = — is the time taken to advect across the edge, 

u 

I 2 

to = is the mean diffusion time for the edge and 

fx = is the time-scale of delivery out of the edge. 
R 

The ratio ^ = ^ is the macroscopic Peclet number for the edge ||39l l60l . If 
I s - 3> 1 then advection is the dominant form of transport across ij, and almost all 
of the material that leaves ij will pass to locations downstream from ij. It is also 
generally true that in the case of high Peclet numbers large concentration gradients 
can persist within each edge ll60l . If rf- <S 1 then diffusion is the dominant form 
of transport across ij, which means that the concentration within ij will tend to 
vary smoothly from node i to j. 

If tT <C tA and tx <C to, then the bulk of the resource will be delivered out of 
the transport network before it transits the edge in question. As a general rule, 
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an efficient transport network will utilize resource over a time-scale tx which is 
similar to the time-scales over which resource transits the whole network. For ex- 
ample, in the case of vascular networks, the oxygen affinity of haemoglobin varies 
with body size, and is related to the circulation time for the species in question 
Il60l 1621 . This makes sense, because if the oxygen affinity of haemoglobin were 
too low for a given body size, red blood cells would become deoxygenated too 
rapidly, and too little oxygen would be carried to the tissues distant from the heart 
and lungs. On the other hand, if a large proportion of the haemoglobin were to 
remain as oxyhaemoglobin throughout the vascular system, only a small fraction 
of the oxygen in red blood cells would be transported to the surrounding tissue. As 
the diffusion coefficient of oxygen is 2 x 10~ 3 mm 2 s _1 ll60l and the velocity of 
flow in a capillary is about 1mm s _1 [ 15], a capillary of length 1 mm has t& = 1 s 
and to = 500 s. Furthermore, as oxygen is delivered throughout an entire network 
of capillaries, it follows that tr 3> tj±. 



2.4 Advection, diffusion and delivery in Laplace space 

As we explain in the Appendix, the quantity of resource in each edge is determined 



by the fundamental Equation (23 ) and the concentrations at the nodes. In Laplace 
space this relationship has a simple algebraic form (Appendix Section 5A). Fur- 
thermore, the Laplace transform of the concentration at the nodes is related to the 
Laplace transform of the net current passing through each node (Section 5B). We 
can invert any solutions that we find in Laplace space back into the time domain 
(Section 5C), and we can also tackle the case of non-zero initial conditions (Sec- 
tion 6). The key idea is that over any time step, the resource in a given edge either 
reaches one of the nodes at either end of the edge, or it remains within the given 
edge. Furthermore, the distribution of resource that has not reached either node is 
equivalent to the distribution of resource that would occur if both nodes were ab- 
sorbing boundaries. These insights enable us to formulate an efficient algorithm for 
calculating how the spatial distribution of resource changes over time in a fixed, ar- 
bitrary network (Section 7). This algorithm couples a network based methodology 
(Sections 7A and 7B) with analytic solutions for individual edges (Section 7C). By 
repeated application of this algorithm, we can also find the spatial distribution of 
resource in a network where the velocities and cross-sectional areas change in a 
piece-wise constant manner over time (Section 7D). In the Appendix we also con- 
sider the more complex case where the cross-sectional areas vary in a continuous 
manner (Section 8). 

The implementation of these algorithms involves the definition of certain constants 
for each edge (Section 5A). In particular, given any Laplace variable s, for each 
edge ij we let 

OHj(s) = Ju^ + ADi^s + Rij). (7) 
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Note that the Laplace variable s represents a rate, and that aij(s) = ctji(s) is pos- 
itive, and dimensionally equivalent to speed. Roughly speaking, ctij{s) represents 
the speed at which resource travels over the time-scale 1/s, with a correction term 
to account for delivery. Since s and Dij are positive and R^ is non-negative, we 
always find that c%(s) > \%Hj . The value of aij(s) depends on mj, and 

U 2 - 

over most time-scales, but for very short time-scales (s S> ^jj 2 - — Rij) almost all 
the movement is due to diffusion, onj 3> Uij and ctij(s) « ^AD{j(s + Rij). 

We let m represent the number of nodes, we let Cj(s) represent the Laplace trans- 
form of the concentration at node i, and (s) is a term that reflects the quantity of 
resource that is initially in edge ij, and which leaves ij by passing through node i 
over the time-scale 1/s (Appendix Sections IIA and IIB). In matrix form we have 



M(s)C(s) =p(s), where 
C(s) = {C 1 (s),C 2 (s),...,C m (s)} T , 

p( S ) = |T 1 ( S ) + ^/3 li ( S ),...,T m ( S ) + ^/3 mj ( S )| , and 



(8) 



Mij(s) 



ik 



v-ik 
2 



Uik( s ) 

2tanh (wfl) 



-Sjjaij(s)e 2Di i 
2si^(^M) 



if i = j, 



otherwise, 



(9) 



We refer to the matrix M(s) as the propagation matrix, and it contains a row and 
column for each node in the given network. Given M(s) and p(s) we can calcu- 
late C(s) using various efficient algorithms, including the stabilized biconjugate 
gradient method (BiCGStab). In most cases this is the most efficient algorithm 
to use, as our matrix M(s) is non-symmetric and sparse [24]. Finding C(s) is 
the most time consuming step of our algorithm, but once we have found C(s) 
for s = In 2/t, 2 In 2/t, . . . , Q In 2/t, we can apply the Gaver-Stehfest algorithm 
Q]|2l[28]|56l|65 ' 66"] to find the concentration at time t for every node in the network 
(further details given in the Appendix). 



2.5 Alternative methods 



As we outlined in Section 2.2 there is a system of equations which govern the 
changing distribution of resource throughout a given network, where the resource 
in question is subject to advection, diffusion and delivery. There are several meth- 
ods that could be applied to solve such a system of equations, in addition to algo- 
rithms described in the Appendix. We could model the movement of resource by 
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taking a particle based approach, where a large number of particles move across 
the network, and the path taken by each particle is determined probabilistically, as 
is the time taken to travel from one node to the next j52l . 

The problem with such particle based approaches is the challenge of avoiding 
under-sampling in the regions of the network that contain a low concentration of 
resource. This problem occurs because, in a finite simulation, the low probability 
paths are, of course, less well sampled, but the fact that such regions are part of the 
network may exert a significant effect on the movement of resource, particularly 
on the higher moments of the transit-times for particles moving across the network 
ll23l l39l l50l . Indeed, that is why the dispersion of tracers can be used to probe 
the structure of networks, and why tracer dispersion plays such a critical role in 
geophysical surveying techniques Il42ll48ll52ll . 

Another possible approach is to employ a finite difference scheme. However, in 
a network where the transport velocities vary over several orders of magnitude, 
straight forward applications of such an approach are not efficient. The problem is 
that the time-scale for updating the concentrations is essentially determined by the 
fastest edge; for stability the distance travelled by advection per time step must be 
smaller than the spatial resolution (ie. the Courant number must be less than one). 
Using such a small time step may be very inefficient in the slower moving regions 
of the network (201 EH. 

3 Vascular geometry and nutrient delivery 

3.1 Calculating the total rate of glucose delivery in idealized vascular 
networks 

We now consider a simple model of glucose moving through a vascular network, 
where the glucose is 'consumed' or transported out of the network by glucose 
transporters on the surface of the vessels. For the sake of simplicity we assume 
that the glucose transporters are uniformly distributed over the interior surface of 
all of the vessels, so the number of transporters per unit length is proportional to 
the radius of the vessel, and the number of transporters per unit volume of blood is 
inversely proportional to the radius of the vessel. 

The rate of glucose delivery reflects the frequency of interaction between glucose 
and the glucose transporters. The kinetics of glucose passing through a transporter 
is rapid ll43l . so high concentrations of glucose are required to saturate the trans- 
porters. Throughout this section we assume that the glucose concentration is below 
the carrying capacity (K m ), and we make the simplifying assumption that the re- 
action rate is proportional to the concentration of glucose and the concentration of 
glucose transporters. In other words, we consider the case where the local delivery 
rate per unit of resource Rij is inversely proportional to the radius of the vessel. 
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We are interested in the total rate of glucose delivery in different networks of cylin- 
drical tubes, as this quantity corresponds to the total rate at which glucose is trans- 
ported out of the vasculature and into the surrounding tissue. We compare different 
network geometries by assuming they have one inlet and one outlet node (nodes 1 
and 2 respectively). We fix the concentration at node 1, inject some volume of fluid 
F per unit time at node 1, and remove an equal volume of fluid at node 2. Given 
the length and radius of each edge, we can calculate the relative conductances, and 
thereby find the medium-current flowing through each edge. This enables us to 
find the velocities Uy as, by definition, the medium-current in each edge is SijUij. 
Given the molecular diffusion coefficient for the resource in question, the disper- 
sion coefficients Dij can be found by Equation (JT|). 

Numerical simulations indicate that the distribution of resource reaches a steady 
state. At steady state, the total rate of resource delivery must equal the current of 
resource entering the network minus the current of resource leaving the network. 



Furthermore, the fundamental advection, diffusion, delivery Equation (23) tells us 
that at steady-state, 

RijQij + Uij - - = 0. (10) 

ox ox A 

It follows that for each edge there must be a pair of constants A and B such that 

u ij+ a ij „ u ij ~ a ij 



q ij (x) = Ae 2D *i +Be 2D ^ , (11) 



where a™ = J u?- + 4Dj 7 i2j 7 -. (12) 



Whatever current of resource and medium we introduce and remove from the given 
network, the steady state distribution of resource must satisfy Equation ( fT0| ). For 
the sake of simplicity we ignore the process of vascular adaptation whereby vessels 
dilate, contract or become apoptotic in response to fluid flow and the associated 
shear wall stress (3l|47l|49j|34l, but as our algorithm(s) can be applied to networks 
with varying cross-sectional areas, we note that such effects could be incorporated 
into a more complex model. 

Our aim is to compare the efficiency of resource delivery for a range of different 
networks, and we do this by calculating the total rate of resource delivery for a rep- 
resentative steady state flow of medium and resource. To find such a representative 
distribution of resource, we suppose that the concentration at node 1 is a fixed con- 
stant k, and that resource leaves the network by flowing from node 2 into a dummy 
edge 2n (see Fig. [3). If we suppose that the concentration at node 2 is c 2 while the 
concentration at node n is 0, Equation ( [TT] ) implies that 

, -S 2n C2e D 2n S 2n C2 
A = -. and B = ; — , 

- a 2n l -°2n l ' 

1 — e D 2« 1 — e D 2™ 
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where I is the length of the dummy edge 2n. Letting I — > oo, we have 



92n(x) = S 2n c 2 e 2 ^n \ (13) 
In this case the flux of resource flowing out of the network at node 2 is 

J2n{x) 



U2nq2n(x) - £>2n ^<?2n (z) 



= F'c 2 , 

x=0 



where Equations ([T]), ( [T2] ) and ( 13 1 tell us that 

pi r, u 2n + Q-2 




(14) 



Given any network of cylindrical tubes with a specified inlet node 1 and outlet node 
2 (see Fig. [3]), and given a molecular diffusion coefficient D m and a local delivery 
rate Rij for each edge, we can find a spatial distribution of resource that reflects 
the network's efficiency as a transport system, and we can calculate the total rate 
of resource delivery in the given case. In particular, it is instructive to calculate the 
total delivery rate at steady state (denoted C to t), which is equal to the total current 
flowing into the network minus the total current flowing out of the network. We 
note that 

C tot = h(t) + h{t) for very large t, (15) 

and we make a fair comparison between different networks by considering the 
following: 

1. In each case, we assume that Fi(t) = F. In other words, at node 1 we inject 
a volume F of fluid per unit time. 

2. We remove an equal volume of fluid from a node 2, so F 2 (t) = —F. 

3. We assume that the flow of fluid is laminar, so the Hagen-Poiseuille equation 

s 2 - 

holds, and the conductance of each edge is proportional to 

4. Given the relative conductances of each edge, and given that a medium- 
current F enters the network at node 1 and leaves the network at node 2, we 
can calculate the velocites Uij EHUD. 

5. All the edges are assumed to be cylindrical and composed of a single ves- 
sel. As we are given the cross-sectional areas Sij we effectively know the 
radius of each edge, as well as and D m , so we can find the dispersion 
coefficients Dij by plugging these values into Equation ([T]). 

6. We suppose that the concentration at node 1 is a fixed constant k at all times. 
This implies that Ci(s) = k/s. 
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7. For each edge ij, including the dummy edge, we suppose that the delivery 
rate per unit of resource Rij is inversely proportional to the radius of the 
vessel. This reflects the assumption that in each vessel there is a fixed density 
of glucose transporters per unit of surface area. 

8. We suppose that the current of resource leaving the network at node 2 is 
completely determined by the concentration at node 2. More specifically, 



we let I 2 (t) = —F'c 2 (t), where F' is given by Equation ( 14 1. Note that the 
value of F' depends on the same cross-sectional area of the dummy edge, 
and in the following section we assume that in each case, S 2n = S12 (see 

Fig- m». 

For the sake of simplicity we assume that each network is initially empty, 
and we calculate the concentrations and total delivery rate for a time point t 
that is sufficiently large for the system to have reached steady state. 



3.2 Analytic solutions to the total rate of glucose delivery in simple 
vascular networks 

We begin by considering glucose delivery by a single vessel (see Fig. [3^), where 
by definition F = F\(t) = S\ 2 i±i 2 . As we are assuming that C\{s) = k/s and the 
network is initially empty, (3ij(s) = for every edge ij (see All), and Equation 



(42 1 tells us that 



M(s) 



It follows that 



so we have 



k/s \_( T!( S ) 
C 2 (s) J { -F'C 2 (s) 

„ , , -M 2 i(s)k 

C 2 (s) = —, =^ r, (16) 

Tx( a ) = -^MnW- M22(s) + F , )■ d7) 



The approximation ctij(s) w a,j = yufj +ADijRij is arbitrarily accurate for 
sufficiently small s, so for very small s we can substitute M^- (0) for My (s). Hence 



Equation (16) tells us that C 2 (s) oc l/s for very small s, and Equation ( 17i tells us 



that Ti(s) oc l/s for very small s. It follows that for sufficiently large t, 

C t ot = h(t) + I 2 (t) = h(t)-F'c 2 (t) 

- <M„,0 ) + «J«). US, 

Note that the terms F' , u and a may vary with the cross-sectional area S: the 
relationship between C tot and the diameter of vessels is plotted in Fig. [4] To 
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replicate the scales of interest in actual vascular networks, we consider edges of 
length I12 = 1mm, and we let D m = 6.7 x 10~ 4 mm 2 s _1 (the molecular diffusion 
coefficient of glucose in water at body temperature). We also let k = 5mmole/litre: 
a typical value for the concentration of glucose in blood. Finally, we assume that 
there is a fixed number of glucose transporters per unit area, so the local delivery 
rate i?i9 = = 4=. The numerical value of the parameter p reflects the 

density of transporters, and their affinity for glucose. To illustrate the biologically 
relevant case, we set p = 0.05mm s _1 . By choosing this value of p we ensure that 
the concentration of glucose drops significantly between the inlet and the outlet 
nodes, but in the networks we consider the concentration does not drop by more 
than an order of magnitude. 



a) b) 3 c) 





■>o ■->• »->o 

1 2 n i 2 n 

Figure 3: Network structure has a critical influence on the total rate of re- 
source delivery. Adding additional vessels may or may not increase the total rate 
of resource delivery, depending on the extent to which the additional edges change 
the time taken to transit the network. We illustrate this effect by considering the 
total rate of delivery in three simple networks, where a current of medium with a 
fixed concentration of resource is introduced at node 1 , and medium and resource 
leaves the network at node 2 by flowing into the dummy edge 2n. 

To produce the continuous curve shown in Fig. [4] we suppose that the medium- 
current F = 0.002mm 3 s~ 1 regardless of the cross-sectional area S\2- In other 
words, we suppose that u\2 = and find the dispersion coefficient D\2 by 
applying Equation ([T}. We have also plotted the case where the pressure drop 
between nodes 1 and 2 is held constant, rather than the medium-current F (see the 
dashed curve in Fig. |4]). The Hagen-Poiseuille equation states that in the case of 
laminar flow, conductance should scale with the square of a vessel's cross-sectional 
area. In other words, in the case of a single edge, maintaining a constant pressure 
drop is equivalent to setting u\2 oc S\2- Consequently, when the pressure drop is 
held constant and the cross-sectional area is very low, very little resource enters the 
network and the total delivery rate is very low. In the case of fixed current, the total 
delivery rate also drops to zero as the cross-sectional area drops to zero, but in that 
case it is because the mean velocity is inversely proportional to the cross-sectional 
area, so when the cross-sectional area is very small, the velocity of flow is very 
large and a relatively large fraction of the resource leaves the vessel without being 
delivered out of the transport network. 

We now consider the other networks illustrated in Fi g. [3| By assumption C±(s) 



k/s, T2(s) = —F' 6*2(5) and Ts(s) = 0. Equation (42i relates these terms to the 



unknowns Ti(s), 62(5) and 6*3(5). As in the previous example, this relationship 
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Figure 4: Total rate of resource delivery in a single vessel. In the case where 
the medium-current is fixed, the velocity of flow is inversely proportional to the 
vessel's cross-sectional area S. Where the pressure drop is fixed, the velocity of 
flow is proportional to S. Our parameters are such that the velocity is lower in 
the case where the pressure drop is fixed, up until the point where S = 1000/um 2 , 
where in either case the velocity u = 2mm s -1 , the time-scale of advection is 
0.5 s and the time-scale of delivery is 0.63 s. Since our parameters imply that the 
medium velocity is smaller in the case of a fixed pressure drop, we also have a 
lower mean concentration and a lower total rate of resource delivery. Note that 
the numerical solution was generated by sampling six points in Laplace space, and 
applying the Gaver-Stehfest algorithm. 

enables us to calculate C tot = h(t) — F'c2{t). In the case of a triangular network 
we let Si 2 = S13 = S32 = 1000/xm 2 and I12 = 1mm (see Fig. [3])). To illustrate 
the effect of short-cuts we vary the length of Z13 = Z32, and see how it effects C to t 
(see Fig. [5]>. As before, is determined by Equation (|T]), k = 5mmole/litre 
and we set Rij = 7== . The velocities Ujj are calculated in two different ways: 

we either fix the total medium-current through the network, or we fix the pressure 
drop between nodes 1 and 2. In either case, the total rate of resource delivery is at 
a maximum when the alternate route is of intermediate length (see Fig. [5]>. 

As a final example we find C to t = h(t) — F' 02(f) in the case where our network 
contains a dead-end (see Fig. ^p). We let S12 = 1000/im 2 , I12 = 1mm, U12 = 
lmms" 1 and U23 = 0, while Dij is determined by Equation ([Tjl. In this case 
we vary the length of the dead-end to see how it effects C tot . As the presence of 
dead-ends vessels can only increase the mean transit time for resource crossing the 
network, we find that increasing the length of the dead-end regions increases the 
total rate of resource delivery (see Fig. [6]). 
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Figure 5: Total rate of resource delivery in a network with alternate routes. 

The medium-current passing through each route will be proportional to its conduc- 
tance (see Fig. [3Jd). If there is a very short alternative route, its conductance will 
be very small, the mean transit time will be very small, and so that the total rate 
of resource delivery will also be small. If the alternative route is sufficiently long, 
most of the resource entering the alternate route is consumed. Further increases in 
the length of the alternate route will decreases the total rate of resource delivery, as 
the medium-current will decrease and so too will the current of resource entering 
the alternate route. It follows that for a fixed current or a fixed pressure drop, the 
total rate of resource delivery is at a maximum for some intermediate length of 
alternate route. 

3.3 Biomedical implications of altering vascular geometry 

Despite their simple nature we now suggest that the results of Section .2 could 
have biomedical implications. Tumours require access to blood vessels for growth 
and metastasis. Consequently anti-angiogenic drugs, which disrupt and inhibit 
the formation of new blood vessels, are a promising avenue for the treatment of 
cancer. As single agents, anti-angiogenic drugs have only produced modest clini- 
cal improvements, but in combination with chemotherapy, the drug bevacizuab (a 
monoclonal antibody against vascular endothelial growth factor) has produced an 
unprecedented increase in survival (5 months) in colorectal cancer patients |[33l . 
This is somewhat paradoxical, as previous studies have indicated that destroying 
the vasculature severely compromises the delivery of oxygen and therapeutics, pro- 
ducing hypoxia that renders chemotherapy and radiotherapy less effective ll35l . 

Tumours are subject to an unusually high interstitial fluid pressure, which may 
collapse blood and lymph vessels, and inhibits the interstitial transport of drugs 
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Figure 6: Total rate of resource delivery in a network with a dead-end. If 

the dead-end region is short (see Fig. |3£), its presence increases the total rate 
of resource delivery by an amount that is proportional to both the volume of the 

dead-end region and the delivery rate per unit of resource -R23. As we assume that 

_ i 

R23 oc S 23 2 , it follows that for sufficiently short dead-end regions the increase 
in the total rate of resource delivery is proportional to I23 y/Syz- The total rate of 
resource delivery reaches a maximum when the time taken to diffuse the length of 
the dead-end region is much greater than the time-scale of delivery. 



ll60l . Furthermore, the blood vessels within tumours are relatively leaky, tortuous, 
and arranged in a haphazard, irregular pattern of interconnection, which results 
in velocities of fluid flow that vary spatially and temporally in a random manner 
(5J [T3l [33 |49j |60l. In healthy tissue the endothelial cells of the vasculature are 
supported by cells known as pericytes, but in tumours the pericytes are loosely 
attached or absent liT71l3"5ll64ll . Anti- angiogenic drugs may impact on drug delivery 
in several ways: they can induce the regression of the particularly leaky vessels that 
lack pericytes [35 ], they can encourage the maturation of the remaining vessels into 
less leaky, less dilated, less tortuous vessels with a greater coverage of pericytes 
lfT7ll3"5ll64l . they can alter the pattern of vascular adaptation (3l|34l|47l and they can 
reduce the interstitial fluid pressure Il34ll54ll60l . Finally, by reducing the number of 
vessels, anti-angiogenic drugs alter the topology of the vasculature. These effects 
have been described as 'vascular normalization' ||5]|35l, and they help to explain 
why the use of anti-angiogenic drugs can actually increase the delivery capacity of 
the vasculature in tumours, increasing the chemo-sensitivity of the tumour itself. 

Our approach helps to illuminate the impact of changes in vascular geometry, as 
we can use our algorithm to compare the delivery rates of various substances for a 
pair of networks (before and after vascular pruning, say). If the delivery rate per 
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unit of resource R is small enough, almost every particle that enters the network 
will exit the network over a time-scale smaller than l/R. In that case the con- 
centration of resource will be approximately constant throughout both networks. 
This implies that any reduction in the total volume of blood vessels will reduce the 
delivery capacity of the network for the substance in question, as the total rate of 
resource delivery is equal to the total volume of blood times the mean concentration 
of resource times R. On the other hand, if R is sufficiently large, almost all the re- 
source that enters the network will be consumed. Again we find that any reduction 
in the volume of blood vessels will reduce the total rate of resource delivery C toU 
but in this case it is because C to t is approximately equal to the current of resource 
entering the network, and reducing the number of blood vessels will increase the 
hydraulic resistance of the network, thereby reducing both the medium-current and 
the current of resource flowing into the network. 

The interesting case is also the most biologically relevant one: where R is such that 
a significant amount of resource is present in the blood that is leaving the network, 
but the concentration of resource entering the network is significantly greater than 
the concentration of resource leaving the network. In this intermediate case, reduc- 
ing the total volume of blood vessels may increase or decrease the delivery capacity 
of the network (that is, the total rate of resource delivery). If we ignore the impact 
of vascular pruning on interstitial pressure, Fig. [6] indicates that removing dead- 
ends can only reduce the delivery capacity of the vascular network. Essentially, 
removing such dead-end regions does not affect the amount of resource entering 
the network, but it does decrease the mean transit time. It follows that the resource 
flowing through the network is more likely to exit the network before it is con- 
sumed, which is to say that removing dead-end regions will decrease the delivery 
capacity of any given network. 

The effect of removing vessels that are an integral part of the network is more com- 
plex. In general, removing the shorter routes between the arteries and veins will 
increase the delivery capacity of the network, as, in the absence of short cuts, any 
resource that enters the network will be forced to spend longer within it, increasing 
the probability that any given particle will be consumed. As an extreme example, 
when an arteriovenous malformation is formed (that is, an abnormal connection 
between arteries and veins) the total volume of blood vessels increases, but such 
a malformation will effectively short-circuit the capillaries in the region, so the 
current in the capillaries and the rate of glucose and oxygen delivery drops dra- 
matically |27l . As Fig. [5] indicates, delivery is optimal when the various routes 
through the vasculature are of similar length, which indicates the importance of 
mechanisms that regulate the demarcation of artery-vein boundaries. This helps to 
explain the importance of Eph/ephrin signals, and other molecular cues that effec- 
tively identify endothelial cells as arterial or venous even before they are fused into 
a functioning circuit H QjO . 

In conclusion, the effect of vascular pruning on glucose delivery will depend on 
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the network structure, and the topological location of the vessels that are pruned. If 
anti-angiogenic drugs eliminate dead-end vessels, the treatment will decrease the 
mean transit time of blood flowing through the tumour. This will tend to reduce 
the total rate of glucose delivery and the chemo-sensitivity of the tumour (though 
this effect may be swamped by other effects of anti-angiogenic drugs, such as a 
reduction in interstitial pressure). On the other hand, if anti-angiogenic treatment 
eliminates the shorter routes by which blood transits through the tumour, our model 
suggests that the effect will be an increase in the total delivery rate of glucose, and 
an increase in the chemo-sensitivity of the tumour. 

4 Concentration in a growing fungal network 
4.1 Modelling the currents in fungal networks 

Multi-cellular organisms need to supply individual cells with the resources neces- 
sary for survival, but while transport in animals and plants is relatively well studied, 
surprisingly little is known about transport in the third major kingdom of multicel- 
lular life. The fungal body or mycelium can be understood as a network of fluid 
filled tubes or hyphae, which grow by osmotically drawing water from their sur- 
roundings while adding material to the cell wall specifically at the tips of the grow- 
ing hyphae Q|45]]. Diffusion may be sufficient to sustain short-range local growth 
when resources are abundant, but foraging fungi such as Phanerochaete velutina 
can grow hundreds of millimeters away from a food source over metabolically inert 
surfaces |[T0l[T6ll46l . Together with various forms of experimental evidence, this 
observation strongly suggests that long-distance transport mechanisms are required 
to deliver nutrients to the growing tips at a sufficient rate, though there are many 
open questions concerning the mechanism(s) of transport lfT2l[T6ll32ll36*ll46l . Vesi- 
cles moved by motor proteins, contractile elements and carefully regulated osmotic 
gradients have all been proposed as mechanisms for driving long range transport 
in fungi lfT6ll36ll46l . Though a fundamental physiological question, which (if any) 
of these mechanisms is important remains debated. 

We note that the fluid within fungal networks is incompressible, and as the network 
grows, there is water uptake in and near the inoculum. It follows that there is a 
mass flow from the sites of water uptake to the sites of growth 11321 . and as the 
tips of the hyphae expand, the cytosol within the organism moves forward along 
with the growing tips [41]. In this section we investigate the argument that this 
form of growth induced mass flow is sufficient to supply the growing tips with the 
resources they require. We do this by modelling advection, diffusion and delivery 
over empirically determined fungal networks. 

To obtain a sequence of digitized fungal networks, we placed a woodblock inocu- 
lated with P. velutina in a microcosm of compacted sand. The growing mycelium 
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was photographed every three days, and the sequence of images was manually 
marked to record the location of nodes or junctions, as well as the presence or 
absence of edges in the fungal network. These edges were not sufficiently well 
resolved to make direct measurements of their diameter from the digitized images. 
However, the reflected intensity, averaged over a small user-defined kernel at ei- 
ther end of the edge, correlated well with microscope-based measurements of edge 
thickness. The observed relationship between image intensity and thickness was 
therefore used to estimate edge thickness across the whole network ifTUll . 

The edges in our fungal networks are composed of bundles of hyphae and trans- 
port vessels bounded by an outer rind ll25ll . Unlike individual hyphae, the edges 
(or cords) in a fungal network have tough hydrophobic coatings which insulate 
them from the environment fT6l l36ll . We make two simplifying assumptions: we 
suppose that all the water and other materials which form the mycelium ultimately 
originate from the inoculum, and we suppose that each edge is composed of trans- 
port vessels, each of which has a typical radius of Qfim [25 ]. Note that the latter of 
these assumptions implies that the hydraulic conductance of each edge is propor- 
tional to its cross-sectional area, as the number of transport vessels in each edge is 
proportional to its cross-sectional area. 

Since the mycelium is composed of incompressible material, the rate of increase 
in the volume of each edge must equal the volumetric rate of flow into that edge 
minus the volumetric rate of flow out of that edge. Together, these assumptions 
enable us to identify a unique medium-current for each edge, namely the set of 
medium-currents that are consistent with the observed changes in edge volume, and 
which also minimize the work required to overcome viscous drag [32]. In effect, 
we simply consider the mycelium as a network of resistors connecting the sources 
of material (the inoculum and shrinking edges) to sinks (the growing edges). This 
enables us to identify a minimal set of growth induced mass flows, and in this 
section we explore whether these currents are sufficient to deliver the resource that 
is required at the growing tips. 

4.2 Modelling resource uptake and delivery 

To find the distribution of resource that results from a given set of currents, we must 
make some assumptions about the rates of resource uptake and delivery. From the 
beginning of each experiment, the inoculum is filled with wood-degrading hyphae, 
so we assume that resource enters the network at the inoculum (node 1) at a con- 
stant rate h(t) = K. The rate of water uptake at the inoculum corresponds to the 
total rate of growth, so our assumptions imply that the amount of water entering the 
network per unit of resource is proportional to the rate of growth. We also suppose 
that throughout the network there is a constant rate of delivery per unit of resource 
R. In other words, where Q(t) denotes the total amount of resource in the network, 
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we suppose that 

j t Q{t) = K - RQ{t). (19) 



As Q(0) = 0, Equation ( 19 1 implies that Q(t) = § (l - e~ Rt ). 

The assumption that K and R are constants implies that the total quantity of re- 
source in the network accumulates over a time-scale , and approaches a steady 
state ^. Furthermore, in our experimental set-up the fungal network attains a max- 
imum volume as there is a finite quantity of resource for the fungi to consume. As 
a final assumption, we suppose that resource accumulates over a time-scale that 
is equal to the time-scale of growth, so that over the course of the experiment the 
mean concentration is approximately constant. More specifically, we let Vp denote 
the maximum volume attained by the mycelium, and we measure the time r that 
elapses before the mycelium attains a volume \Vp. We then assume that Q(t) is 
half the maximum quantity of resource. The numerical value of K reflects the units 
we use to measure the concentration, so without loss of generality we can assume 
that the mean concentration at time r is 1. It follows that Q(r) = ^ and ^ = Vp. 
This implies that 

log(2) Vplog(2) 
R=^^- and K= - SV ' , (20) 
r r 

so we have Q(t) = Vp(l — 2~). 



4.3 Modelling the spatial distribution of resource in empirical net- 
works 



To apply our minimal model for the distribution of resource in a growing fungal 
network, we require empirical values for Vp (the maximum volume attained by the 
network) and r (the time taken to grow to volume \Vp). We also require the ad- 
jacency matrix of the network, the lengths lij and the cross-sectional areas Sij(t n ) 
for each edge ij and each time point ii, . . . ,ijy 

For each time interval, the first step is to calculate the unique set of medium- 
currents which are consistent with the observed changes in volume, and which 
minimize the work required to overcome viscous drag [32 ]. We suppose that over 
the time interval t n < t < t n+ \ the cross-sectional area Sij(t) = \{Sij(t n ) + 
Sij(t n+ i)). Furthermore, as the edges are composed of a bundle of transport ves- 
sels, we suppose that the conductance of each edge is proportional to its cross- 
sectional area. Finally, we calculate whether each of the nodes is a source or a 
sink. Where Fi denotes the net medium-current flowing out of node i, we let 

{— Ylj^i Fj ^ n °de i i s me inoculum, 

(21) 
^ sM-s tJ (t +1 ) otherwise _ 
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If the edges around node i are growing, then Fj is negative and node i is a sink, 
which is to say that more medium-current flows into the links of node i than flow 
out. If the edges around node i are shrinking, or if node i is the inoculum, then 
Fi is positive and node i is a source. Circuit theory tells us that we can use the 
conductance of each edge and the net current flowing out of each node to determine 
the pressure difference between any pair of nodes iPTl l32l . Furthermore, given 
the conductance of edge ij and the pressure drop between nodes i and j, we can 
uniquely determine the medium-current Fij(t) for each edge in the network. This 
medium-current is constant over the time interval t n < t < t n+ \, and it does not 
depend on the constant of proportionality between the cross-sectional area of the 
edges and the conductance of the edges. 

The edges or cords in a fungal network have a complex structure [25 1 , and mass 
flows occur in transport vessels that occupy some fraction A of the cross-sectional 
area of each edge. The medium-current in an edge is equal to the mean velocity of 
flow times the total cross-sectional area of the transport vessels, so for each edge 
and each time interval we have 

Ui . (t) = 2i ^' (t) . (22) 

X[Sij(t n ) + Sij(t n+ \)j 

As we wish to investigate whether growth-induced mass flows are sufficient to 
carry resource from the inoculum to the tips over the time-scale of growth, we set 
A = 1. In other words, given values for the medium-currents, we let the veloci- 
ties of mass flow be as small as possible by maximizing A. Also note that, given 
the observed changes in volume, and given the assumption that resource and water 
only enters the network at the inoculum, the medium-currents that we identify are 
as small as possible (in the sense that any other set of medium-currents consistent 
with the observed growth would require more work to overcome viscous drag). We 
are interested in finding the distribution of a generic source of energy and carbon, 
so we let D m = 6.7 x 10~ 4 mm 2 s -1 (the molecular diffusion coefficient of glu- 
cose), as this is representative of the diffusion coefficient of a small molecule. We 
also assume that in each edge the advection and diffusion of resource occurs within 
some number of transport vessels of radius 6/mi ll25ll . so once we have found the 
mean velocity of flow Uij(t), we can use Equation (jl]) to find the piece-wise con- 
stant dispersion coefficient Dij (t) . 

The delivery rate per unit of resource is assumed to be the same in every edge, and 



the value of Rij = R is given by Equation ( 20 1. In each experiment the parameters 
Vf and r were chosen to ensure that, over each time step, the mean concentration 
is as close to 1 as possible. The nutrient and water content of woodblocks can 
vary, resulting in more or less vigorous growth. In the first replicate we found that 
Vf = 393mm 3 (20% of the volume of the woodblock) and r = 242 hours. In the 
second we found that Vp = 372mm 3 (19% of the volume of the woodblock) and 
r = 468 hours. In the third we found that Vp = 616mm 3 (31% of the volume of the 
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woodblock) and r = 367 hours. Finally, the Laplace transform of the net quantity 
of resource leaving the inoculum is assumed to be Ti(s) = — — , and for every 
other node Tj(s) = 0. We now have all the parameters we need to implement 
either of the algorithms described in the Appendix. That is to say, we can calculate 
the spatial distribution of resource that would arise if the cross-sectional areas of 
the edges vary in either a step-wise or continuous manner, where the volumetric 
currents are determined by the measured changes in volume. Whichever algorithm 
we we employ, at each time t n we record the resulting spatial distributions of re- 
source by dividing each edge into Nij line segments such that < 1mm. These 
quantities of resource per unit length are then treated as an initial condition over 
the following time step, and the concentrations at time t n are identified by dividing 
qij(x,t n ) by Sij(t n ). 



4.4 Results of the simulation 

Three fungal networks were grown and digitized, and the observed changes in 
fungal volume were used to determine the minimal currents consistent with the 
changing volume, as well as the uptake rate and decay rate of a generic form of 



resource (see Sections |4. 1 [|4.2| and |4.3| >. In each of three experiments our model 
suggests that the growth induced mass flows were sufficiently large to spread the 
resource from the inoculum out to the growing tips over the time-scale of growth 
(see Figs. [7]and[8]>. 

This result is somewhat counter-intuitive, as in most of the edges the mean velocity 
of the growth induced mass flows is very low |32|. Indeed, if we pool the data from 
all three experiments and over all time steps, 70% of the edges have a mean velocity 
that is so small that over the course of one week, resource travelling at that velocity 
would move less than the 20mm that resource would typically travel by diffusing 
in one dimension (75% of edges in Experiment 1, 59% in Experiment 2 and 64% in 
Experiment 3). Over the time-scale of two hours, only 4% of edges have a velocity 
great enough to carry resource further than the 2.2mm that is typically travelled by 
diffusion alone (2% in Experiment 1, 6% in Experiment 2 and 6% in Experiment 
3). 

Despite the modest scale of the advection in most of the edges, the fraction of edges 
in which the mean velocity is significant suffices to spread the resource from the 
inoculum out to the growing tips (see Figs. [7]and[8]). We also calculated the distri- 
bution of resource that results if the cross-sectional areas Sij (t) vary continuously 
over each time step (see AIV), but the results were almost identical to the simpler 
case where the cross-sections are varied in a stepwise manner. 
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Figure 7: Concentration of resource in Experiment 1. Experiment 1 after 6 days 
(a, d, g and j), 12 days (b, e, h and k) and 18 days (c, f, i and 1). Diagrams (a - c) 
illustrate the concentration of resource that would occur in the absence of advec- 
tion, where resource enters at the inoculum at a rate of K = 1.125^mole hour -1 , 
and r = 242 hours. Diagrams (d - f) illustrate the concentration of resource where 
fluid and resource enter the network at the inoculum, and the medium-currents are 
consistent with the observed changes in volume, while minimising the work re- 
quired to overcome viscous drag (see Section 4.1 1. As before, resource enters at 
the inoculum at a rate of K = 1.125/umole hour -1 , and r = 242 hours. Note that 
at any point in time, the concentration near the tips can be greater than the concen- 
tration near the inoculum. This is possible because resource enters the network at 
a constant rate, but the rate of water influx at the inoculum corresponds to the total 
volumetric rate of growth. Consequently, as the total volumetric rate of growth 
increases, the concentration of resource in the fluid near the inoculum decreases. 
In (d), for example, the fluid in the tips contains more resource than the fluid near 
the inoculum, but when that fluid first entered the network (at the inoculum) it con- 
tained an even higher concentration of resource. We cannot directly measure the 
delivery rate of resource, so to assess the sensitivity of our model to the parameter 
R, we also consider the cases where we half and double the delivery rate R. Dia- 
grams (g - i) illustrate the concentrations that occur when the medium-currents and 
rate of uptake are as before, but the local delivery rate has been halved. Diagrams 
(j - 1) illustrate the concentrations that occur when the medium-currents and rate of 
uptake are as before, but where the local delivery rate has been doubled. 
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Figure 8: Concentration in the tips relative to the concentration elsewhere. In 

each of the three experiments, the concentration at the tips (nodes of order one) 
was larger than the mean concentration in the network as a whole. 
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Figure 9: Network structure and the efficacy of growth induced mass flows as 
a means of transport. Two contrasting examples of networks with growth induced 
mass flows. 

4.5 Discussion of growth induced mass flows 

In a fungal network, the incompressibility of aqueous fluids ensures that growth 
in one part of the network requires the presence of fluid flows in the supporting 
mycelium. By controlling the spatial location of growth, maintaining the appropri- 
ate turgor pressure and by forming cords that are insulated from the environment, 
fungi can ensure that there is a long range flow of fluid from the sites of water up- 
take to the sites of growth ll32l . Furthermore, the structure of the network is critical 
for ensuring that growth induced mass flows can carry resource from the inoculum 
to the growing tips over a reasonable period of time. In this regard, it is instructive 
to compare a growing linear network to a growing branching tree (see Fig. [9]>. 

Suppose that the tips in both of the model networks illustrated in Fig. [9] grow a 
unit distance from the inoculum per unit time, and that each edge has unit length 
and volume. In the case of a linear network, and in the absence of diffusion, it 
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will take n units of time for the resource to travel from the inoculum to the nth 
edge. It follows that if the time-scale of delivery is n, then growth induced mass 
flows in a linear network cannot supply resource over length scales greater than 
n. In the case of a branching tree the volume of the nth generation is greater than 
the total volume of all the preceding generations (see Fig. [9]>, so it will take less 
than a unit of time for resource to travel from the inoculum to the nth generation. 
Provided that the concentration of resource at the inoculum remains sufficiently 
high, there is no limit to the size of branching tree that can be filled with resource 
by growth induced mass flows, even if the local rate of resource delivery is high. 
Also note that in the absence of diffusion, the fluid exactly at the growing tips 
is never replaced by the fluid entering at the inoculum. This implies that growth 
induced mass flows alone cannot supply resource to the growing tips: diffusion 
and specific transport mechanisms are essential for transporting resource across 
the newest generation of edges. 

Returning to our model of transport and resource delivery in an empirical fungal 
network, it could be argued that the total rate of resource delivery in an edge should 
be proportional to the volume of that edge, rather than being proportional to the 
quantity of resource contained within that edge. However, as our model results in a 
fairly constant concentration throughout the network, changing from first order to 
zero order delivery is unlikely to make a significance difference to the concentra- 
tion at the tips (unless, of course, we change the mean amount of time that elapses 
between resource entering the network and the resource being consumed). It might 
also be argued that the growing hyphal tips are responsible for a significant fraction 
of the resource consumption |[T2ll46ll . After all, as material is added to the growing 
cell wall, the concentration of that material in the cytoplasm must be depleted near 
the region of growth. Although our model does not include a term for consump- 
tion due to growth, it does indicate that growth induced mass flows are sufficient 
to carry resource across the network over the time-scale of growth. Furthermore, 



Equation ( 22 1 indicates that if the growth induced mass flows are confined to trans- 
port vessels that only occupy a fraction A of the total cross section of each edge, 
then the mean advective velocities will be greater than our minimal estimates by a 
factor of 1/A. 

While advective mass flows carry resource over long distances from the inoculum 
out towards the growing tips lfl6l l36l 1461 . diffusion and active transport mecha- 
nisms may be essential near the sites where the cell wall is expanding. This follows 
because the cytosol within the hyphae moves forward at the same rate as the grow- 
ing tips fill , but to transport resource from the base of the hyphae to the growing 
tips, the resource has to move faster than the rate of growth. Complex cellular 
machinery regulates the addition of material to the cell walls, ensuring that the 
growing hyphae exhibit polar growth, and only expand at the hyphal tips f71l45l. 
We note, however, that our model suggests that vesicles carried by motor proteins 
or other active transport mechanisms may not be needed for longer range transport 
within fungal networks. 
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Conclusion 



In the Appendix we present two algorithms for calculating the concentration of 
resource that arises when a given material is subject to advection, diffusion and 
local delivery out of the transport network. The resulting spatial distribution will 
depend on the time-scale required to transit the network, and the time-scale of de- 
livery (see Section [23] >. Nature is full of networks in which materials within a fluid 
are transported by advection and diffusion while being consumed or delivered, so 
these algorithms have many potential applications. In particular, our modelling 
framework can be applied to the case of glucose delivery through a vascular net- 
work. By analyzing simple, idealized vascular networks we found that the total 
rate of glucose delivery depends on the network structure (see Section [3]), and in 
some cases increasing the volume of blood and the number of glucose transporters 
can actually decrease the total rate of glucose delivery. This counter-intuitive re- 
sult can occur because the additional vessels can decrease the time taken to transit 
the network, allowing a greater fraction of the glucose to pass through the network 
without encountering a transporter. Finally, we employed our algorithms to im- 
plement a model of transport in a growing fungal network (see Section [4]). The 
expansion of fluid filled vessels requires the movement of fluid, and in three empir- 
ically determined fungal networks we found that the minimum currents consistent 
with the observed growth would effectively transport resource from the inoculum 
to the growing tips over the time-scale of growth. This suggests that the active 
transport mechanisms observed in the growing tips of fungal networks may not be 
required for long range transport. 
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Appendix: Mathematical methods. 



5 Advection, diffusion and delivery in Laplace space 



As we explain in the Main Text, we are motivated to consider the case where re- 
source is lost or delivered out of a given network at some local rate, while the 
resource that remains within the network moves by advection and diffusion. Such 
a process will result in a spatial distribution of resource that changes over time. We 
only consider longitudinal coordinates along each edge ij, using real numbers x to 
denote distances from node i, where < x < kj. Each edge contains a quantity 
of resource, which must satisfy the one-dimensional advection-diffusion-delivery 
equation 

dt + ijqij + '' dx ij dx 2 



o. 



(23) 



where qij is the quantity of resource per unit length, mj is the mean velocity, Dij 
is the dispersion coefficient and Rij is the rate at which a unit of resource is lost, 
or delivered out of the network. As we explain in the Main Text Section 2.2, Fick's 
law implies that 



m = 



Uij(t)qij(x,t) - A 



dqjj(x,t) 
' dx 



(24) 



s=0 



where denotes the net quantity of resource leaving node i at time t. Note that 
by the conservation of mass any resource that enters node i must leave node i (as 
nodes have zero volume), so Ii(t) = unless resource is entering the network at 
node i. Nodes where resource enters the network are referred to as inlet nodes. 

Given such a system of fundamental equations, we want to find the quantity of re- 
source throughout the network in an efficient manner. It is convenient to follow the 
approach of [39], which exploits the Laplace transform. This operation effectively 
weights the different time-scales over which resource may move from one node to 
another, and it is an efficient way of coping with the wide range of velocities that 
our network may contain. In particular, we take advantage of the following prop- 
erties of the Laplace transform C(qij(x, t)) = qij(x, t)e~ st dt = Qij(x, s): 



£( d gtJ Q f ' t ) = sQij(x, s) - qij(x, 0), and 



dt 



-£««<...>■ 



If Uij and Dij are constant over time, it follows from Equation ( 23 1 that 



(s + Rij)Qij + 



dQ^ 
dx 



D„ 



d 2 Qj 
' dx 2 



(25) 
(26) 

(27) 
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Furthermore, Equations (24 1, p5| ) and ( [26] ) imply that 
where Tj(s) denotes the Laplace transform of Ii(t). 



d 



x=0 



(28) 



5.1 Zero initial conditions in an edge 



We begin by considering an initially empty edge, before extending our results to 
the more complex case of non-zero initial conditions. We let qij(x, 0) = 0, and 



consider the homogenous case of Equation (21): 



(S -|- Rij)Qij ~t~ Uij 



&Qij ® Qij 

dx %3 dx 2 



0. 



(29) 



By solving this ODE in the usual manner, we find that for some pair of constants 
A and B, 



Qij(x } s) — Ae 
where otij(s) 



u ij+ a ii( s ) 



+ Be 



ufj + ADi^s + Rij). 



(30) 
(31) 



Note that the Laplace variable s represents a rate, and that a>ij(s) = ctji(s) is pos- 
itive, and dimensionally equivalent to speed. Roughly speaking, «y(s) represents 
the speed at which resource travels over the time scale l/s, with a correction term 
to account for delivery. Since s and Dij are positive and R,^ is non-negative, we 

always find that aij(s) > uy . When s = j^ 2 - 



Rij , Equation (31) implies that 



a: 



2uf, ; . The value of ctij(s) depends on u^, Dy and R^ over most time 



scales, but for very short time scales (s 3> ^g 2 - — Rij) almost all the movement is 
due to diffusion, aij 3> and ay (s) « \/4Dy(s + 



Equation ( |30| ) tells us that for any positive number s, we can find ^4 and i? and 
express Qij(x, s) in terms of the quantity of resource at either end of the edge. For 
any given s, we denote the quantity of resource at the ends of each edge by 



Xijis) 



Qij(0,s) 



and 



Xji(s) = Qji(0,s) — Qij(kj,s). 



(32) 



For each edge ij, it is convenient to define two, dimensionless ratios between time 
scales: 



Uijlij 

2Du 



and hij(s) 



2D 



(33) 



1.1 



Setting x = and x = kj tells us that 



Xij = A + B 



and 



Xji = Ae^ 9ij+hi ^ + B e (f«-M. 



(34) 



29 



A, B, Xij, Xji, ctij and hy are all functions of the Laplace variable s, but this 
dependence is omitted for the sake of clarity. Equation ( 34 1 tells us that 

^ e (9y+M + (Xij - A)e^~ h ^ = X, 
Xjie-Ba - Xjj< 



so 



A 



2 sinh(/ij 



and likewise 



B 



X ije h v - X,;< ■< • 



2 sinh(/tj 



(35) 
(36) 



Note that if Uij is negative, the medium-current flows towards node i and the 
macroscopic Peclet number for the edge 

ij is = -2gij [39, 60]. Assum- 

ing that edge ij is initially empty, we can find Qij(x, s) by substituting Equations 
([33]>, ([35]> and ([36]) into Equation ([30]>, giving us 



X; 



sinh(' 



sinh(/i 



-e *■» +a 



sinh( I |^ J ) 
&wh(hij) 



e '« 9lJ . (37) 



5.2 Advection, diffusion and delivery in an initially empty, static net- 
work 



Having examined the case of a single edge, we now turn to the problem of coupling 
the edges of a network such that the concentrations vary continuously as we move 
from one edge to another. For each node i we have Cj(s) = J °° Ci(t)e~ st dt. 
Assuming that the cross-sectional areas Sij are constant, Equations (Main Text 3), 



and ( 32 1 imply that for all edges ij we have 



Ci(s) 



Xjj(s) 

Si-j 



and Cj(s) 



Xjj(s) 



(38) 



Enforcing this equation ensures that the Laplace transform of the concentration at 
node i is consistent for all edges ij, ik, and so on. In general, we may not know the 
Laplace transform of the node concentrations C{s) = {Ci(s), . . . , C m (s)}, where 
m is the number of nodes. However, given I(t) = {h(t), . . . , I m (i)} (the net cur- 
rent of resource leaving each node), we can calculate T(s) = (Ti(s), . . . , Y m (s)} 
(the Laplace transform of I), and, in the following manner, calculate C(s). If we 
substitute Equation ([30]) into Equation ([28]), noting that x = tells us that 



T<(s) 



Y j ( f{B-A) + V f{A + B). 



Equations ( 33 1 and ( 34 1 imply that A + B = Xij , and 
1 



B — A 



2 sinh(hij) 



Xije h v - .\j,r " ■ - .V ; , ( ''• + X ije ~ h ^ 



30 



so we have 



2 sinh(/ijj^ 



cosh(/ij 



+ 



(39) 



Equations ( 38 1 and ( 39 1 imply that 



Ti(s) 



E 



C{(s)Si. 



u 



+ 



a; 



2 2tanh(/tj 7 ) 



C,( S )S, 



2 sinh(/ijj 



• (40) 



In other words, for each node i we have a linear equation in Ci(s), C2(s), . 
Hence where C(s) and T(s) are column vectors, we thus have 



where 



M( S )C( S ) = T(a), 

Hk 



My(s) 



+ 



Oik 



2tanh(/i ifc ) 



2 sinh 



{his) 



if i = j, 



otherwise. 



i C m (s). 
(41) 

(42) 



We refer to the matrix M(s) as the propagation matrix, and it contains a row and 
column for each node in the given network. Given M(s) and T(s) we can calcu- 
late C(s) using various efficient algorithms, including the stabilized biconjugate 
gradient method (BiCGStab). In most cases this is the most efficient algorithm to 
use, as our matrix M(s) is non-symmetric and sparse [24]. 

Equation pT) implies that the diagonal elements M(s) are all positive. Further- 
more, Mjj(s) = if and only if there is no edge between i and j, and the other 
off-diagonal elements are negative. We note that if there is resource at node j, it 
may be transported along ij, bringing resource to i and reducing Tj(s) (the Laplace 
transform of the net current flowing out of node i). Resource can only reach node i 
along the edges ij, ik, etc, so Tj(s) is completely determined by the concentration 
at i and the concentrations that flow through the nodes adjacent to i. As Tj(s) 
is the Laplace transform of the net current flowing out of node i, and resource at 
nodes j ^ i can flow into node i, the off-diagonal elements of M(s) are negative, 
and zero if i and j are not directly connected. 

Multiplying |My-(s)| by Cj(s) gives us the Laplace transform of the current of 
resource flowing from j to i, so roughly speaking, Mjj(s) represents the size 
of the volumetric current from j to i, over the time scale 1/s. Note that if Uij is 
positive, then the medium-current flows from i to j, Mjj(s) < Mjj(s) , and 
there is a greater flow from i to j than the other way around. That is to say, when 
the medium-current is from i to j, the value of Cj(s) has a greater influence on 
the value of Yj(s) than the influence of Cj(s) on the value of Tj(s). Also note 
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that the ratio of My (s) to Mjj(s) depends on the Peclet number "jj ' j = 2^, as 



Mjj(s) : Mjj(s) is equal to 1 : e 



> Uij 



For very short time scales we have a very large s, and by Equation (31 
and a.ij ^/ADijS. In this case the off-diagonal elements of M are very small, 
and Mjj « Sik^rf- ~ X^fc Sik\fDikS- In other words, over very short time 
scales resource is lost from the nodes by a process of diffusion, but it does not 
have time to reach the other nodes. Over longer time scales the difference between 
and ctij is smaller, the off-diagonal elements of M are larger, and effect of 
advection is greater. 



5.3 Inverting from Laplace space 

We now have a method for finding the Laplace transform of various quantities, 
and in this section we consider how to transform these quantities into the time 
domain. More specifically, we have seen that for a given Laplace value s, we 



can find M(s) and T(s). We can therefore use Equation (41 1 to find C(s) = 
{Ci(s), . . . , C m (s)}, the Laplace transform of the concentrations at each node. 
Furthermore, we can use Equation ( |37| ) to calculate Qij (x, s) in terms of the bound- 
ary conditions Xij(s) = £(^(0, tJV. and Xji(s) = C[qij(lij,t)). In other words, 
for each edge and each Laplace variable s, we can find an algebraic expression for 
Qij(x, s) in terms of the boundary conditions Xij(s) and Xji(s), but we have yet 
to show how we can numerically invert such quantities into the time domain. 

As we can calculate any sequence of real valued sample points in Laplace space 
and we wish to calculate the corresponding value at a given point in time, it is ap- 
propriate and efficient to apply the Gaver-Stehfest algorithm [2, 22]. The key idea 
behind this algorithm (and other, related algorithms) is the notion of construct- 
ing a sequence of linear combinations of exponential functions, in order to form a 
weighted delta convergent sequence [2, 22, 65, 66]. That is to say, we consider a 
sequence of functions 5 n (x, t) such that for any function q that is continuous at t, 
we have 

poo 

5 n (v,t)q(v)dv = tq n (t), (43) 



o 



where q n {t) — > q(t) as n — > oo. As we shall see, there are weighted delta conver- 
gent sequences of functions such that 5 n (v , t) is of the form 

n 

Sn(v,t) = ^uj i e-r L , (44) 
i=i 

where > for all i, and the terms 6{ and Ui do not depend on t. Now, if we 
suppose that our function q does not increase exponentially, then the Laplace trans- 
form Q{s) = J °° e~ sv q{v)dv is well defined for all positive numbers s. Hence the 
existence of Q(s) for all positive s is a reasonable assumption, given the context 
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in which our functions q arise. Assuming that Q(s) is well denned for all positive 



numbers s, Equations (43 1 and (44i imply that 



Qn(t) = - S^uJie^ q{v)dv = - VwiQ(^). 

1 J ° i=i f i=i 1 

Gaver [28] employed the sequence of functions 

„ , * , (2n)\ , v In 2 , „ , v In 2 . „ 

S n (v,t) = In 2 1 ' 1 - e-— * e~— n , 
n!(n — lj! 

but the resulting terms g n (i) converge to logarithmically slowly. Gaver also 
showed that the quantity q n (t) — q(t) can be expanded in terms of inverse powers 
of n, which enabled him to accelerate the convergence of his original sequence of 
approximations [28]. The most useful formula for finding an accurate estimate of 
q(t) based on a linear combination of the Gaver estimates was derived by Stehfest 
[56], who stated that 

In 2 In 2 

l{t) ~ Quit) = — 2^ K nQ{n— ), where (45) 



=i 



K = , ir+ n/ 2 mm( f /2) 

n y ' ^ (n/2 -k)lk\(n-k)\(2k-n)V 

fe=[(n+l)/2] y ' ' v ; v ; 

and is even. Note that the terms n n can be extremely large, and that the value of 
K n depends on the parameter Q. Furthermore, increasing the parameter Q increases 
the accuracy of our estimate q(t) w q~n(t), provided that we have sufficient system 
precision to utilize the exact values for n n . 

The Gaver-Stehfest algorithm is very efficient and accurate, but it requires high 
system precision for the weights n n if it is to yield accurate estimates for q(t). In- 
deed, if we wish to produce an estimate of q(t) that is accurate to N significant 
digits, we must calculate the values of n n with an accuracy of about 2.5iV signifi- 
cant digits [1, 2]. Fortunately, to calculate q(t) accurately we do not require such a 
disproportionately high level of accuracy in the values of Q(s). 

If the transform Q(s) has all its singularities on the negative real axis, and if the 
function q(t) is infinitely differentiable for all t > 0, extensive experimentation [1, 
2] indicates that the relative error 

q{t) - qn(t) 



10 -o-45n (46) 



provided that the values n n have been calculated with sufficient precision [1, 2]. 
If the function q does not satisfy the above conditions qa(t) may converge to q(t) 
rather more slowly, but as a rule of thumb setting Q, = 10 and using standard 
double precision for the weights K n will ensure that the Gaver-Stehfest algorithm 
produces inversions that are accurate to at least three significant digits. 
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6 Non-zero initial conditions 



6.1 Non-zero initial conditions in a single edge 



We now consider advection, diffusion and delivery along a single edge ij, where 
the initial condition q%j{x, 0) is non-zero. We let the length of ij equal /, the lon- 
gitudinal dispersion coefficient is D, the local delivery rate is R and the mean 
velocity is u. 

We have seen that for any positive Laplace constant s, Qi(x, s) = e^ 9+h ^ and 
Q2(x, s) = e^ 9 ~ h ^ satisfy the homogeneous analog, Equation (29). Furthermore, 
the Wronskian 



ur i \ /-> t s.dQ%{x,s) dQi(x,s) 
Wij[x,s) = Qi(x,s) — — Q 2 (x,s) 



dx 



dx 



D 



By the method of variation of parameters, 



f[x,3,qij(y,0)) 



a 



o 



a; 







is a particular solution to the fundamental Equation (27 1. 



(47) 



Note that f(0,s,q) = for all initial conditions q. Also note if q = q\ + q 2 
then/(x,s,g) = f(x, s, qi) + f(x, s, q 2 ). Since f(x, s, qij(y, 0)) is a particular 
solution of Equation ( [27] ), for each edge ij there is a pair of constants A and B 
such that 



Qyfas) =Ae^+ h n+ Be^- h n+ f( x , s, qij (y, 0)). 



(48) 



Because /(0, s, q) = for all initial conditions q, Equations (33 1 and (48 1 imply 
that 



X ij = Q ij (0,s) = A + B, 



and 



(49) 



Xji^Qijik^s) = Ae^+ h )+Be^ + f{l,s,q tj (y,0)). (50) 



We can therefore express A and B in terms of Xij and Xji. Indeed, substituting 
Equation (49) into Equation (|50]> and multiplying both sides by e~ 9 tells us that 



X jie - 9 = A(e h -e- h )+X ije - h + e- 9 f(l,s, qil (y,0)). (51) 



We let 



-ae 



2 sinh(/i) 



f{l,s,qij(y,0)), 



(52) 
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and explain its physical significance in Section 6.2 
Equations ([49]) and (BTT) imply that 



A 



B 



2 sinh(/i) 
2 sinh(/i) 



+ 



ft 



and 



a 



fa 
a 



(53) 



and substituting Equation (53 ) into Equation (48 1 tells us that for any initial condi- 
tion qij(y,0), 



Qij "SJ 



Xj{€ ^ X{j€ 

2 sinh(/i) 



a 



(54) 



f X ije h - Xjje-9 A 
\ 2 sinh(fo) a 



+ f{x,s,qij(y,0)). 



we 



6.2 Non-zero initial conditions over a network 

Having analyzed the case of a single edge with a non-zero initial condition, >. 
now consider an entire network, and find an exact solution that ensures that for all 
t > 0, the concentration varies continuously as we move from one edge to another. 
The first step in finding this solution is to note that Equation ( 47 ) implies that 



dx 



(u — a) 



2Da 



where for the sake of clarity we drop the subscript ij from Uij, ctij, kj, gij, hij and 

Dij. Note that for any initial condition qij(y, 0), we have df (*' s «f^) \ x=o = o. 
It follows that 

dQij(x,s) 



dx 



i=0 



/% u + a f X jie ~9 - X ije - h 
D 2D \ 2sinh(ft,) 



+ - 



u — a ( Xije h — Xjie i 



2D V 2sinh(/i) 



(55) 



Now, recall that Tj(s) denotes the Laplace transform of the net current of resource 
flowing away from node i, and that Tj(s) = unless i is an inlet node. Substituting 



Equation (55 ) into Equation (28 1 gives us 



2 2tanh(/ti ? ) 
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Assuming that the cross-sectional areas Sy are constant, Equations (38) and (56) 
imply that 



T 4 (a) = C i (s)J2S ij (f + 



(X; 



2 2tmh(hij) 



(57) 



In matrix form we have 

M(s)C(s) =p(s), where (58) 
C(s) = {C 1 (s),C 2 (s),...,C m (s)} T , 

Pi (s) = T t ( S ) + J2^(s) (59) 



and M(s) is the propagation matrix as in Equation (42 ). Note that the effect of the 
initial conditions on the concentration at the nodes is completely captured by the 
terms and that, as before, the propagation matrix M(s) relates the concen- 

trations at the nodes to the net currents flowing out of the nodes. Furthermore, by 



comparison with Equation (42), we can see that the concentration at the nodes is 
the same as would be the case if the network were initially empty, and the Laplace 
transform of the net current leaving node i were Pi(s) rather than Tj(s). 

In effect, the formalism of the propagation matrix enables us to substitute an initial 
condition in the edges around node i for a boundary condition at node i. For 
each node i and each Laplace variable s, this boundary condition is of the form 
Ylj Pij( s )- Intuitively speaking, the term /%(s) represents the Laplace transform 
of the quantity of resource that first leaves edge ij by arriving at node i. Note that 
we have not calculated the impact of the initial condition q%j(x, 0) on the future 
concentration profile qij(x, t) for t > 0: we have simply calculated the impact of 
the initial conditions on the concentrations at the nodes (see Section [7]). 

Since ctij(s) 3> Uij and hij(s) 3> gij for large s, for very short time steps t we 
have sinh(fey) 3> max [e 9ij , e~ 9ij ~\ . It follows that over short time scales, the off- 
diagonal elements of M(s) will be very small. If the entries in the ith column of 
M(s) are very small, it may be numerically difficult to calculate Cj(s), as any error 
in our estimate for Q(s) would have very little impact on the value of M(s)C v (s). 

In practice this is not a significant problem, as when we solve the above system of 
linear equations to identify we make the initial guess that Cj(s) = Cj(0)/s, 

which would be the correct value if the concentration at node i remained constant. 
For numerical reasons we may not be able to identify the exact value of Cj(s), but 
this problem only arises when the bulk of resource around node i does not leave the 



edges around node i over the time scale 1/s. As we shall see in Section 7.2 those 
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are precisely the circumstances under which the value of Cj(s) has little impact on 
our calculation of the spatial distribution of resource within the edges at a given 
time t ps 1/s. 

7 Efficient calculation of resource distribution 

If we wish to find the concentration at various points in the network other than the 
nodes, there are two ways we can proceed. The first method is to treat each point of 
interest as an additional node. The problem with this approach is that it increases 
the size of the propagation matrix, and finding C(s) by inverting the matrices M(s) 
is the major computational cost of the propagation matrix algorithm. Furthermore, 
although this approach can be used to find the exact concentration at each of a given 
set of points, it does not provide a means of finding the exact quantity of resource 
between a given pair of points. We could approximate the total quantity of resource 
between two points by assuming that the concentration varies in a linear manner 
from one point to the next, but as the exact solution may contain boundary layer 
effects, we might require a very high spatial resolution to ensure that such a linear 
approximation is accurate. 

A different approach, which we take, provides an exact solution for the total quan- 
tity of resource within each section of the network, regardless of the spatial reso- 
lution. The key conceptual step involves partitioning the resource into two parts. 
Strictly speaking our approach is mathematically continuous, but we can imagine 
that the resource is composed of particles, which either leave or do not leave a 
given edge over a given time scale. We let q%j{x, t) denote the quantity of resource 
per unit length at the point < x < kj in edge ij and time t, where a given 
particle only contributes to qij(x,t) if it has passed through a node (any node) 
by time t after initialization. More precisely, we work in Laplace space and let 
£{<iij( x it)) = Qij( x , s )- This term denotes the Laplace transformed concentra- 
tion profile that would occur if the network was initially empty, and if the Laplace 
transform of the net current leaving each node was Pi(s) = Tj(s) + J^j Pij( s )> 
rather than Tj(s). 

As we have seen, the impact of the initial condition on the concentration at the 
nodes is completely captured by the constants Pij(s). However, Qij(x,s) and 
qij(x, t) do not fully capture the influence of the initial condition qij(x, 0) on the 
concentration profile qij(x,t) for t > 0. In addition to qij(x,t) (the quantity of 
resource that has reached a node over the time scale t), we must also consider the 
resource that starts in edge ij, and which does not reach node i or j over the time 
scale t. We let qij (x, t) denote the quantity of such resource at the point < x < kj 
in edge ij and time t, where by definition 

qij(x,t) = qij(x,t) - qij(x,t). (60) 
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We can calculate the concentration at each node by calculating Pij(s) for every i 
and j, and by using the propagation matrix. Furthermore, because at time none 



of the resource in edge ij has had time to reach a node, we can apply Equation 37 
and find Qij(x, s) in terms of the boundary conditions Xij(s) and Xji(s). Given 
Qij(x, s) for s = ^p, . . . n lj j^, we can apply the Gaver-Stehfest algorithm and find 
qij (x, t). In addition, we solve a separate PDE for each edge, which tells us how the 
resource that stays within each edge has evolved over a given time step t. That is to 
say, for each edge ij we find qij(x, t), given that qij{x, t) satisfies the fundamental 



advection-diffusion-delivery Equation (23 1, qij(x, 0) = % 0), qij(0, t) = and 



qij(kj,t) = 0. Finally, Equation 60 tells us that qij(x, t) = qij(x, t) + <jij(x, t). 



In particular, we consider the case where the initial condition is stepwise constant, 
and edge ij is divided into Nij sections of equal length. We let k^j (t) denote the 
mean quantity of resource per unit length in the nth section at the given time t, 
where by convention the first section is next to node i and the iVfh section is next 
to node j. For any t > 0, we can employ the following algorithm to find an exact 
solution for the updated mean quantities per unit length, 

k%\t) = ^l I ' ' ,/,,!,-. /),/,-. (61) 



7.1 Stepwise constant initial conditions 

We are interested in calculating how the quantity of resource in a network changes 
over time, given that the resource decays and is subject to advection and diffusion. 
In particular, it is convenient to consider a stepwise constant initial condition, as 
we can then calculate how the total quantity of resource in each segment of the net- 
work has changed by time t. The first step in this calculation is to find the Laplace 
transform of the concentrations at each node C(s). As we have seen, to calculate 
C(s) we must first find Mjj(s) and T(s), which do not depend on the initial con- 
dition. For each sample point s and each edge ij we must also calculate Pij{s) and 
Pji(s), which capture the effect of the initial condition qij(x, 0). In particular, we 
start this subsection by considering the case where the initial condition is 



qij(x,0) 



h If n-l i. . < T x 
otherwise, 



where n < N. We will find our solutions for other initial conditions by summing 
the solutions for various initial conditions of this form. For the sake of clarity we 



drop the subscripts ij from lij, Nij, gtj and hij, and note that Equation (47) tells 
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us that for this initial condition 

f(l,s,qij) = 



ke g+h r%l 



a; 



<9+hV Jdx 



+ 



a: 



^fdx, 



OLijiUij + Oij) 



(62) 



Substituting into Equation ( |52| ) yields 

ft; 00 



ke N n 9 



4(s + Rij) sinh(/t) 



N — n t j h_ 

e n n ( eN - e n ){ aij - Uij 



)(< 



(63) 



Recall that f(x,s,qi + 92) = f(x,s,q\) + /(x, s,^)- Since Equation (52 1 is 
linear, it follows that if the initial condition contains several blocks of resource, 
each block makes its own separate contribution to /%(s) and (3ji(s). Let xq = 
0, x\ = jj, X2 = ■ ■ ■ , xn = I, and suppose that for all 1 < n < N we have 



qij(x, 0) = A;-™* 1 for all x n _i < x < x n . 



(64) 



Given such a stepwise constant initial condition, we can calculate /%(s) by sum- 
ming the contribution of each of the blocks of resource. That is to say, Equation 



( [64] ) becomes 



N 



J 4(s + sinh(/t 



+ e " "y( e w -e f ) (o^- + n^j 



(65) 



We can find /3ji(s) by using the above formula, substituting — ^ for gji, —u^ for 
u,i and n+l ^ for fcI-7 . It follows that 

= ^4(/+ J R lj )sinh(/ iii ) 



"~ n h ■ I -21- !H\ / \ 

e at ft 'j( e f - e Jv j ( aij + Uij ) 



+ e f 



e n 



e 



)(« 



(66) 
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7.2 Resource that leaves its initial edge 



If a particle leaves edge ij and reaches node % or j over the relevant time scale, 
it contributes to /%(s) or f3ji(s), and hence it contributes to our solution Cj(s), 
Cj(s) and C(qij(x,t)) = Qij(x,s). On the other hand, at time none of the 
resource has reached the nodes, so the initial condition q~ij(x, 0) = 0. It follows 



that if the cross-sectional areas are held constant, we can apply Equation (37 1. In 
other words, we can find Qij(x,s) by effectively considering an initially empty 
network, where resource is introduced at the nodes at a rate which exactly matches 
the rate at which resource reaches the nodes in the case where the network has the 



given non-zero initial condition. Equation ( 37 1 also accounts for the impact of any 



inlet nodes, in the case where resource is being added to the network. 



We can therefore use Equations (58 1, (65 1 and (66 1 to find C(s) = {Ci(s), . . . , C m (s 
and in the case where the cross-sectional areas are constant, we can use Equations 
(37 1 and (38) to express Qij(x,s) in terms of the boundary conditions Xij — 



SijCi(s) and Xp = SijCj(s). Since £( f qij(x,t)dx) = f Qij(x, s)dx, we can 



find J qij(x,t)dx by calculating f Qij(x, s)dx for s 
applying the Gaver-Stehfest algorithm. 



\n2/t, . . ,,N1n2/t and 



We suppose that edge ij is divided into Nij sections of equal length, and for the 



sake of clarity we drop the subscripts ij from Dij, Uj and Nij. 



We let y^\t) 



denote the mean value of qij(x, t) in the nth section of edge ij, and note that by 
definition 



N 

T 



(jij(x, t)dx. 



(67) 



Defining Y^ l \s) = C[y^\t)) we have 



n r¥ 



Qij(x, s)dx 



NX, 



'j 



I sinh(/i,; 



I — x N 



smh(hij — : — )e 9ij i dx 

1 7 I 



NX 



I sinh(/ijj) jn-i 



sinh(/ijj -) e giJ * i dx 
l <■ 



ND 



I sinh(/i 



Xije^-X 3i e **_ (gij .- feij)f 



u 



a;. 



+ 



i.r 



ua + a 
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which implies that 



Yij{s) = Vij(s){aij + Uij) x 

. n-JV re + iV ^ n-JV-1 n + JV-1 fe , 

( g jv jv w — g i\r a« jv "»j ) 

+ //,;(•-)(<>,, - Uy) X 

/ 2JV-n+l i ti 2N—n r . 

Xjjfe AT 9l J JV "*J — gAfS^ JV "'Jj 



L (9ij+^ij) 



Nije h v 



where = — ^ 7 — r- 

4(s + i2»j)Zjj smh{hij) 



(68) 
(69) 



7.3 Resource that remains in its initial edge 

Over the time scale t, not all of the resource will leave the edge in which it started. 
To find qij(x, t), the quantity of resource that has not left edge ij, we must solve 
the advection, diffusion, delivery problem for each separate edge ij, where nodes 
% and j are absorbing boundaries and the initial condition qij(x, 0) = qij(x,0). 
The resulting solution accounts for those particles which do not reach a node in the 
relevant time-scale. In particular, we consider the case where the initial condition 



is stepwise constant, as in Equation (64i. 



The fundamental Equation ( 23 1 tells us that for each edge 



0_ 
dt 



d 2 







^ij q x 2 ?*j Ui Q x Qij Rij Qij ■ 



(70) 



Furthermore, we are looking for a real valued function such that qij(0,t) = and 
q~ij{kj,t) = for all t. These conditions imply that we can express q~ij(x, t) in the 
following form: 



qij{x,t) 



where A"? 



e e lJ sm (— — ). 



m=l 



— ~\~ R-ii 

4A, 13 



(71) 



The parameters A m can be found by taking Fourier transforms. More specifically, 
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we know that q~ij(x, 0) = % 0), so 



J2 A m si 



.mirxs 



sin 



n=l 



J sin(— — ) sin(- — Jax 

Jo Hj Hj 



qij(x,0)e J! « and 
r if m 7^ n, 



-y- if m = n. 



It follows that for every positive integer m, 



sm( — — Jgij(x,0)e « ax. 



In particular, consider the case where the initial condition is stepwise constant, and 
of the form described by Equation (64). We have 



A 1 



i, m V fr (n) 

Vi] K ij 

71=1 



,mirx^ 



(I I I I \ iXj \ / I 

— — j — cos(- 

ij ' Hj 



JLI. . 



where 



u" 1 
Pi 3 



n=l 



e 'J k 



m ( u (n+i) _ k (n) \ ( 9ij „ inf mmr 



sin 



* J\irm v N, 



+ cos 



'J 



Mil 



(72) 



(73) 



We are now in a position to find 
as Equation (J7T]) implies that 



n i 

AT. . /" W- • K? 



''V 



TI — 1 I 

TV- ■ l4 J 



qij(x,t)dx, 



jy. . rw^kj hi* °° 



( ; y A m e x % 

n-l , ^— ' 



sin 



Ida; 



m=l 



m=l 



7rm 



sm 



run-, -^f . ,m(n— l)ir, 



e N{ i sin ( 



AO 



,m(n— l)7r N ,mnir, 
+ ( e <j cos( — J — cos(^ J 



4 j 



(74) 
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Note that a™- — > — as m — > oo, and likewise A m £ 0(m 1 ). In contrast e 
tends to zero much more rapidly. Indeed, we note that 



E 

rn—n> 



\ m i 



< 



( u ii 
e Ud h 



r 

\ 4£>y 



m=n> 



~J2 x 



xe l » dx 



I 2 - 



2Q,'-K 2 D ij t 



(75) 



It follows that the relative error 

>oo A'"t v^£V A"H 



m=l e -Em=l e 



v^oo A^t 



< 



Z^m=n> e 



v^oo XTit 



< e 



whenever we have 



e « < e- 



*J m=l 



E 



A m t 

e ! ^ . 



(76) 



We can therefore be confident that if we truncate the sum in Equation d74b at m 



Q', the relative errors in our estimates for z^\t) will be smaller than e provided 



that Q! satisfies Equation (76 1. Also note that Equation (71 1 tells us that if Dijt > 



9 \ m t / 9 

If a then e *» decreases rapidly, so SI does not need to be large unless D^t <C 
Furthermore, if uf^t 2 > I 2 - then most of the resource will leave edge ij over the 
time scale t, and qij(x, t) will only make a small contribution to the total value of 
qij (x, t). In that case using a small value of fi' will produce very accurate estimates 

fack$\t) even if £> -i < Z£. 



7.4 Calculating the total quantity of resource in each segment of a 
network 

Suppose that we wish to calculate the mean concentration per unit length in each 
segment of a network at time t, such that each part of our final answer has a relative 
error e < io~°- 45r2 , where Q, is an even integer. The first step is to set s = Q In 2/t, 



and apply Equations (65 ) and (66 1 to find /%(s) and (3ji(s) for each edge ij. We 
then compute M(s) and and employ the BiCGStab algorithm to find C(sq), 
starting with the initial guess that for each i, 

c '<»"» - = s^jw (77) 
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This initial guess for the value of C (sq) would be correct if the concentration at the 
nodes was constant, and making such a guess can help to speed up the process of 
finding the true value of C(sq). At each step, when we have identified C(s) such 
that M(s)(7(s) = p(s), we store the vector C(s) and repeat for s = sn_i, . . . , s±, 
where s n = n In 2/t. The only difference is that for subsequent applications of the 
BiCGStab algorithm, we can take advantage of the approximation 

77 + 1 

Ci( Sn ) « -—0^+!). (78) 

n 

This is generally a better initial guess than that provided by Equation ( [77] ), so the 
BiCGStab algorithm converges on the solution more rapidly. Given Ci(s n ) and 



Cj(s n ), we can use Equation (68 1 to calculate Y^J n \s n ) for each section in the 
edge ij. Having found Y^ n \s n ) for each 1 < n < f2, we can apply the Gaver- 
Stehfest algorithm to obtain y\j(t), and we repeat this process for each edge in 



the network. Finally, for each edge ij we can use Equations ( |7lj ), ( |73| ) and ( |72| ) 



\ rn f 

to calculate a sequence of values for e « , jjff^ and A m until we reach an integer 



il such that e » satisfies Equation (76 1. We then employ Equation (74i to find 



z!-j\t), . . . , {t) (the mean quantity of resource in ij that has not reached a 
node), and note that for each section of the network the mean quantity of resource 
per unit length 

k$\t)= V W(t)+z$\t). (79) 

Unless there are many sections in each edge, finding the vectors C(s) such that 
M(s)C(s) = p(s) is the most time consuming step of the computation, as it ef- 
fectively involves inverting an m x m matrix M(s), where m is the number of 



nodes. We also note that Equation (42) implies that if hij is larger than 10 (say), 
then the matrix M(s) may be close to singular, making it computationally difficult 
to calculate C(s). Fortunately this problem is easy to avoid, as we can simply in- 
troduce an additional node k at the midpoint of edge ij. This increases the size of 
the matrix M(s), but the lengths 1^ and l^j will be half the length l^. As we have 
seen, the ratio Mjj(s) : Mjj(s) is equal to 1 : e 29ij , so adding additional nodes 
greatly reduces the ratio between the entries of M(s), and can make it significantly 
easier to find the vector C(s). 

Finally, we note that this algorithm can be adapted for the case where the cross- 
sectional areas Sij(t) vary continuously over time (see Section IV). However, even 
in the case where Sij(t) varies continuously over time, our method requires that 
over each time step the lengths lij, mean velocities mj, decay rates Rij and dis- 
persion coefficients D{j are held constant. In the case where we wish to find the 
concentration of resource in a changing network, we simply vary all the param- 
eters in a stepwise manner, finding the spatial distribution of resource at the end 
of each time step, and treating that distribution as an initial condition for the fol- 
lowing time step. In the case of the fungal networks that we analyze in the Main 
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Text, this approach yields very similar results to the more complex algorithm with 
continuously varying Sij (t) that we outline in the following section. 

8 Advection, diffusion and delivery in a changing network. 

We now consider the case where each cross-sectional area Sij (t) varies monotoni- 
cally over time, while the lengths kj, mean velocities Uij, decay rates Rij and dis- 
persion coefficients Dij remain constant. Equation (1) tells us that the dispersion 
coefficients Dij will be constant and equal to the molecular diffusion coefficient 
D m if the edges are sufficiently narrow or the velocities sufficiently low. Alterna- 
tively, D,^ and would remain constant despite the changing cross-sections in 
the specific, but biologically relevant case where the edges are composed of some 
variable number of tubes of fixed radius r^, while the pressure at each node re- 
mains constant over time [26, 27, 44]. In that case the conductance of ij will be 
proportional to Sij(t), so the medium-current at time t will be proportional to the 
pressure drop times Sij(t), and the velocity will be constant. 

More generally, if we are considering advection, diffusion and delivery over a net- 
work where the parameters u^, Sij, Rij and change over time, it is reasonable 
to assume that u^, Rij and are piece-wise constant provided that the time 
scales for transiting the edges ij are small compared to the time scales over which 
u^, R^ and are changing. For example, in the case of vascular networks the 
cross-sectional areas of capillaries, the velocity of blood flow and the rates of re- 
source delivery may vary over time, but such changes occur over time scales that 
are large compared to the time it takes to transit a capillary. In such a case Rij 
might represent the local rate of glucose delivery per unit of glucose in the blood 
(for example), and the following algorithm enables us to calculate the concentra- 
tions that arise at times t\,t2, etc, as we vary Sij(t) in a continuous manner, while 
, R^ and Dy vary in a stepwise manner, being held constant between each of 
the time points of interest. 

Now, suppose that we want to know how the spatial distribution of resource in a 
network changes over a time-scale r. We let SV,(t) denote the area of edge ij 
at time t, and where <%(0) and <%(r) are given quantities, it is mathematically 
convenient to set 

(t) = {2Sij (r) - (0)) + (2Sij (0) - 2S l3 (r)) e^ 1 

= aij + bije^ 1 . (80) 

By adopting this functional form for Sij (t), we are assuming that the cross-sectional 
areas Sij vary in an approximately linear manner over the time scale of interest r. 
In fact, the rate of change 4 Sij halves over the time scale < t < r. Also note 
that if the given cross-sectional areas S^ (0) and Sij(r ) are non-negative, it follows 
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that Sij(t) > for all < t < r. However, as t oo, Sij(t) 2%(r) - <%(0), 
which may be negative. Since the value of Sij(t) has little effect on our calcula- 
tions for t > t, we are not introducing a major source of error when we allow the 
possibility of negative values for 5y (t) at time points beyond the time of interest. 

8.1 Propagation matrices for a changing network 

Suppose that the cross-sectional areas <SV/(t) are of the form described by Equation 



( 80 1, and Uij and Dij are constant. By definition we have (x,t) = (t)cij (x,t), 
so taking Laplace transforms gives us 

Qij(x,s) = C(aijCij(x,t)) + C(bije r t c i j(x,t)) 



= aijCij(x,s) +bijCij(x,s + ln2/r). 
In particular, writing s' = s + In 2/t, we define 

Xij(s) = Qij(0, s) = aijCi(s) + bijCi(s') and 
Xji(s) = Qij(kj, s) = aijCj(s) + bijCj(s'). 

Substituting Equations ([82]) and ([83]) into Equation (39) tells us that 



a 



2 2tanh(/i ii ) 



Yl (oijCj(s) + bijCj(s'] 



ctije 9 ^ 



2 sinh(hy) ^- 

In matrix form we have the equivalent of Equation (41): 

V(s)C(s) + W(s)C{s + In 2/t) = p(s), 
where = {d( S ), C 2 (s), C m (s)} J , 

Pi(s) = Ti(s) + E Pij(s), 



(81) 

(82) 
(83) 



(84) 



(85) 



and 



W -(s) 



2 2tanh(h lfc ) 



2 sinh 



u ik I "it 

2 2tanh(/j ifc ) 



2 sinh 



if i = j, 
otherwise, 

if i = j, 
otherwise. 



(86) 



(87) 
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8.2 Calculating C(s) in a changing network 



To find the boundary conditions Xij(s) and Xji(s) for each edge, we must first 
calculate the Laplace transform of the concentration at each node. Furthermore, 
in order to apply the Gaver-Stehfest algorithm, we need to calculate C(s) for s = 
si, . . . , sq where s n = nln2/r. As in the case where the cross-sectional areas 
of the network remain constant, we can calculate V(s), W(s) and p(s) for any 



positive integer s. Since s n +i = s n + In 2/r, we can use Equation (85 1 to relate 
the vectors C{s n +\) and C(s n ). That is to say, for each n we have 

V(s n )C(s n ) = p(s n ) - W(s n )C{s n+1 ). (88) 

To begin this iterative process of finding C{s n ) from C{s n+ \), we must first es- 
timate the value of C(sq) for some integer Q. In finding such an approximation 
the first point to note is that the value of Ci(sq) is predominantly determined by 
the value of Cj (i) over the time-scale < t < ^ : the value of a (t) for t > ^ is 
relatively inconsequential. 



The second point to note is that by Equation (80), the cross-sectional area at time 

SijQ = ^•(0)(2% i -l)+^(r)(2-2^). 

If Q = 13 (say), we have ~ 0.9,% (0) + 0.1%(r). We can estimate 

Ci(sn) by assuming that the cross-sectional area of each edge does not vary over 
time, but remains constant at %(n)- This enables us to apply Equation (41), and 
thereby obtain an estimate for Cj(sn). More specifically, for each node i we make 
the initial guess that 

(an approximation that would hold exactly if the concentration at the nodes re- 
mained constant). We can then employ the BiCGStab algorithm to home in on 
a more accurate solution to M(sq)C , (sq) = p(sa). Once we have obtained an 



estimate for C(sn), we can employ Equation (88 1 to find C(sq-i), . . . ,C(s±). 
Also note that when the relative change in cross-sectional area is small, b{j is small 
compared to etjj, so the elements in the vector W(s n +i)C(s n ) are small compared 
to the corresponding elements in ~V(s n )C(s n ) and p(s n ). This means that when 
the relative change in cross-sectional area is small, any errors in our estimate for 
C(sq) have little effect on the calculated values for C(sq-i). 

8.3 Resource distribution in a changing network 

We now combine the preceding observations, and present an algorithm for calcu- 
lating the concentration in each section of a network with changing cross-sectional 
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areas after a given time r. More specifically, we suppose that for each edge in 
our network the mean velocity ttjj and the dispersion coefficient Dy remain con- 
stant, while the cross-sectional areas Sy(t) vary smoothly from Sij(0) to Sij(r) 



in accordance with Equation (80 1. The first step in our algorithm is to pick an odd 
integer S7. This determines the scale of the errors in our final answer, and the ratio 
of the errors to the true values which will be of the order e = io -0 - 45 ^^ 1 ). As a 
rule of thumb setting 17 = 13 and using standard double precision for the weights 
K n will ensure that the our final answers are accurate to at least three significant 
digits. 

We let sq = rH r n2 and assume that the cross-sectional area of edge ij is held 
constant at 5V/(^). These cross-sectional areas can be used to find a propagation 
matrix M(sn), as described by Equation (25). Furthermore, the given initial con- 
dition can be used to calculate p(sq), as described by Equations (42), (48) and 
(49). 

Equation (41) tells us that M(sq)C(sq) = p(sq), so we are now in a position to 
find C(sq). More specifically, we can employ the BiCGStab algorithm, starting 
with the initial guess that for each i, 

V k {1) 



171n2 v ' nin2X)j5ij(0)" 

If the concentration at the nodes remained constant, this estimate would equal the 
exact solution. Once we have run the BiCGStab algorithm, we store the vector 



C(s n ), and use Equations (42), (48), (49), (j86|> and (|87j» to find V(s„), W(s n ) 
and p(s n ) where n = 17 — 1. 



These matrices are related to one another by Equation (88 1, so once again we can 
employ the BiCGStab algorithm to find C(s n ). As in the case of a network with 
constant cross-sectional areas, the matrices tend to be close to singular if any of 
the terms li ^^( s ) we large (greater than 10, say). This makes it computationally 
difficult to calculate C, but this problem is easy to avoid as we can simply intro- 
duce additional nodes along the edge ij, thereby reducing the lengths iy. Having 
described the network in terms of nodes and sufficiently short edges, we start the 
BiCGStab algorithm with the initial guess that for each i, 

n + 1 

Ci(s n ) ~ Cj(s n .fi). 

n 

This process is iterated until we have found C(sq), . . . ,C(si). Now, for each 



edge and each integer n = 17 — 1, . . . , 1, we can use Equations (82) and (83) 
to find the boundary conditions Xij{s n ) and Xji(s n ). As in the case of constant 
cross-sectional areas, the Laplace transform of resource that has reached a node 
is denoted by Qfj{x, s n ), and this quantity is related to Xy(s n ) and Xji(s n ) by 

Equation (20). We can therefore use Equation (51) to find \s) = C{y$\t)), 
the mean value of QfAx, s) in the n'th section of the edge ij. 
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Having found Y^ 1 ' (s) for each section in the network, we can apply the Gaver- 

Stehfest algorithm to obtain y^\t). Finally, we note that the quantity of resource 
within each edge only changes because of advection, diffusion and delivery. Vary- 
ing the cross-sectional area of an edge will affect the quantity of resource that enters 
that edge, but it will not directly affect the quantity or distribution of resource that 
remains in ij without reaching nodes i or j over the time-scale t. For example, if 
an edge ij is shrinking and the fluid within ij is leaving that edge, it is assumed 
that the effect of this mass flow is entirely captured by calculating the appropriate 
velocity term Uij. Given that this is so, for each edge ij we can use Equations (54) 
and (56) to calculate a sequence of values for e X ^ 1 , fjy and A m until we reach 
. \ n ' t 

an integer U such that e «j satisfies Equation (58). This ensures that when we 
employ Equation (56) to find z^\t), . . . , z^' 3 \t) (the mean quantity of resource 
in ij that has not reached a node), the errors are very small. Finally, we note that 
for each section of the network the mean quantity of resource per unit length 

8.4 Results 

We have presented two algorithms for calculating the concentration in a network as 
it changes over time, with resource subject to advection, diffusion and local deliv- 
ery out of the network. In the first, cross-sectional areas and the medium-currents 
in each edge change in a stepwise manner, while in the second the cross-sectional 
areas and medium-currents in each edge vary continuously. The former algorithm 
is about twice as fast, and experimentation indicates that it is more numerically 
stable. 

By design, the total medium-current to pass through each edge in a given time step 
is the same whether the cross-sectional areas vary continuously or in a stepwise 
manner. Consequently, if the concentration at each node does not change dramati- 
cally over a given time step, the amount of resource to pass through each edge over 
the time step in question will be similar whichever of the two algorithms we apply. 
However, the two algorithms may give very different results for the concentration 
profile within a given edge, if there is a difference in the initial concentration at the 
nodes at either end. 
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Figure 10: Concentration in a linear network and branching tree. In each 
case resource is added at the source at a constant rate, the local delivery rate 
throughout each network is 0.02hour _1 , and the diffusion coefficient is D m = 
6.7 x 10~ 4 mm 2 s _1 . Every hour, both networks grow 1mm further from the source, 
and the branching tree bifurcates every 10mm. The cross-sectional area of each 
new edge either varies in a stepwise manner from to 6/mi (red lines), or the 
cross-sectional area varies continuously from to 6/im (black lines). Note that the 
same amount of resource is contained within the linear network and the branching 
tree, but the branching tree has a greater volume. After 70 hours the mean concen- 
tration in the linear network is over 1 8 times greater than the mean concentration 
in the branching tree, but the concentration near the tips is only about 5 times as 
great. 
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